我如何select纬度和经度在地球仪上形成一个"rectangle"？
How do I select latitudes and longitudes to form a "rectangle" on a globe?
我有一个数据集，其中包含经纬度网格上的值。我需要从这个数据集中 select 绘制北美近乎完美的 "rectangle"。某物 like this，但位于北美上空：
1。如何选择经纬度？
由于经度向两极汇聚，我需要更多向北的经度。
这是我的拙劣且可能不正确的尝试。我猜这是 cartopy 中的单行代码，但我不知道我在寻找什么转换。
我的矩形的高度在北纬 0° 到 75° 之间。我正在计算每个纬度的经度跨度，使得水平宽度（以米为单位）与 0° 纬度时从 215° 到 305° 经度的距离相同。 （矩形水平居中 260°。）
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
def long_meters_at_lat(lat):
"""Calculate distance (in meters) between longitudes at a given latitude."""
a = 6378137.0
b = 6356752.3142
e_sqr = a**2 / b**2 -1
lat = lat * 2 * np.pi / 360
return np.pi * a * np.cos(lat) / (180 * np.power(1 - e_sqr * np.square(np.sin(lat)), .5))
min_lat, max_lat = 0, 75
min_lon, max_lon = 215, 305 # Desired longitude spread at at min_lat
central_lon = (max_lon + min_lon) // 2
dist_betn_lats = 111000 # In meters. Roughly constant
lat_range, lon_range = np.arange(max_lat, min_lat-1, -1), np.arange(min_lon, max_lon+1)
x_idxs, y_idxs = np.meshgrid(lon_range, lat_range)
y_meters = (y_idxs - min_lat) * dist_betn_lats
y_lats = y_idxs + min_lat
dist_betn_lons_at_min_lat = long_meters_at_lat(lat_range[-1])
x_meters = (x_idxs - central_lon) * dist_betn_lons_at_min_lat # Plus/minus around central longitude
x_lons = central_lon + np.round(x_meters/long_meters_at_lat(lat_range)[:, None]).astype('uint16')
assert ((x_lons[:, -1] - x_lons[:, 0]) <= 360).all(), 'The area is wrapping around on itself'
x_lons = np.where(x_lons >= 360, x_lons-360, x_lons)
这就是y_lats, x_lons的样子，看起来很正常（右上角的低经度绕了360°）。
(array([[75, 75, 75, ..., 75, 75, 75],
[74, 74, 74, ..., 74, 74, 74],
[73, 73, 73, ..., 73, 73, 73],
...,
[ 2, 2, 2, ..., 2, 2, 2],
[ 1, 1, 1, ..., 1, 1, 1],
[ 0, 0, 0, ..., 0, 0, 0]]),
array([[ 87, 91, 94, ..., 66, 69, 73],
[ 97, 101, 104, ..., 56, 59, 63],
[107, 110, 113, ..., 47, 50, 53],
...,
[215, 216, 217, ..., 303, 304, 305],
[215, 216, 217, ..., 303, 304, 305],
[215, 216, 217, ..., 303, 304, 305]], dtype=uint16))
2。我如何在地球上绘制这些 latitudes/longitudes 处的数据？
我在下面尝试了明显的，但只是在右边有一条窄条。
crs = ccrs.PlateCarree()
u = np.random.rand(*x_lons.shape)
v = np.random.rand(*x_lons.shape)
ax = plt.axes(projection=ccrs.Orthographic())
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax.set_global()
ax.scatter(y_lats, x_lons, u, v, transform=crs)
plt.show()
明显的错误是代码中 (long, lat) 的反转。
这是要尝试的正确代码。
# (second part only)
crs = ccrs.PlateCarree()
u = np.random.rand(*x_lons.shape)
v = np.random.rand(*x_lons.shape)
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-80, central_latitude=30))
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax.set_global()
ax.scatter(x_lons, y_lats, u, v, transform=crs)
plt.show()
编辑 1
这里是完整的代码，只在地图上的某个矩形形状内绘制数据。
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
import numpy as np
import matplotlib.patches as mpatches
# part 1 (minor change)
def long_meters_at_lat(lat):
"""Calculate distance (in meters) between longitudes at a given latitude."""
a = 6378137.0
b = 6356752.3142
e_sqr = a**2 / b**2 -1
lat = lat * 2 * np.pi / 360
return np.pi * a * np.cos(lat) / (180 * np.power(1 - e_sqr * np.square(np.sin(lat)), .5))
min_lat, max_lat = 0, 75
min_lon, max_lon = 215, 305 # Desired longitude spread at at min_lat
central_lon = (max_lon + min_lon) // 2
mean_lat = (max_lat + min_lat) // 2
dist_betn_lats = 111000 # In meters. Roughly constant
lat_range, lon_range = np.arange(max_lat, min_lat-1, -1), np.arange(min_lon, max_lon+1)
x_idxs, y_idxs = np.meshgrid(lon_range, lat_range)
y_meters = (y_idxs - min_lat) * dist_betn_lats
y_lats = y_idxs + min_lat
dist_betn_lons_at_min_lat = long_meters_at_lat(lat_range[-1])
x_meters = (x_idxs - central_lon) * dist_betn_lons_at_min_lat # Plus/minus around central longitude
x_lons = central_lon + np.round(x_meters/long_meters_at_lat(lat_range)[:, None]).astype('uint16')
assert ((x_lons[:, -1] - x_lons[:, 0]) <= 360).all(), 'The area is wrapping around on itself'
x_lons = np.where(x_lons >= 360, x_lons-360, x_lons)
# part 2
from_lonlat_degrees = ccrs.PlateCarree()
# map projection to use
proj1 = ccrs.Orthographic(central_longitude=central_lon, central_latitude=mean_lat)
u = np.random.rand(*x_lons.shape) # 0-1 values
v = np.random.rand(*x_lons.shape)
# auxillary axis for building a function (lonlat2gridxy)
axp = plt.axes( projection = proj1 )
axp.set_visible(False)
# this function does coord transformation
def lonlat2gridxy(axp, lon, lat):
return axp.projection.transform_point(lon, lat, ccrs.PlateCarree())
fig = plt.figure(figsize = (12, 16)) # set size as need
ax = plt.axes(projection=proj1)
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
# create rectangle for masking (adjust to one's need)
# here, lower-left corner is (-130, 15) in degrees
rex = mpatches.Rectangle( ax.projection.transform_point(-130, 15, ccrs.PlateCarree()), \
6500000, 4500000, \
facecolor="none")
ax.add_artist(rex)
bb = rex.get_bbox() # has .contains() for use later
# plot only lines (x,y), (u,v) if their
# (x,y) fall within the rectangle 'rex'
sc = 1. # scale for the vector sizes
for xi,yi,ui,vi in zip(x_lons, y_lats, u, v):
for xii,yii,uii,vii in zip(xi,yi,ui,vi):
xj, yj = lonlat2gridxy(axp, xii, yii)
# check only p1:(xj, yj), can also check p2:(xii+uii*sc, yii+vii*sc)
# if it is inside the rectangle, plot line(p1,p2) in red
if bb.contains(xj, yj):
ax.plot((xii, xii+uii*sc), \
(yii, yii+vii*sc), \
'r-,', \
transform=from_lonlat_degrees) #plot 2 point line
pass
# remove axp that occupies some figure area
axp.remove()
# without set_global, only rectangle part is plotted
ax.set_global() # plot full globe
plt.show()
我有一个数据集，其中包含经纬度网格上的值。我需要从这个数据集中 select 绘制北美近乎完美的 "rectangle"。某物 like this，但位于北美上空：
1。如何选择经纬度？
由于经度向两极汇聚，我需要更多向北的经度。
这是我的拙劣且可能不正确的尝试。我猜这是 cartopy 中的单行代码，但我不知道我在寻找什么转换。
我的矩形的高度在北纬 0° 到 75° 之间。我正在计算每个纬度的经度跨度，使得水平宽度（以米为单位）与 0° 纬度时从 215° 到 305° 经度的距离相同。 （矩形水平居中 260°。）
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
def long_meters_at_lat(lat):
"""Calculate distance (in meters) between longitudes at a given latitude."""
a = 6378137.0
b = 6356752.3142
e_sqr = a**2 / b**2 -1
lat = lat * 2 * np.pi / 360
return np.pi * a * np.cos(lat) / (180 * np.power(1 - e_sqr * np.square(np.sin(lat)), .5))
min_lat, max_lat = 0, 75
min_lon, max_lon = 215, 305 # Desired longitude spread at at min_lat
central_lon = (max_lon + min_lon) // 2
dist_betn_lats = 111000 # In meters. Roughly constant
lat_range, lon_range = np.arange(max_lat, min_lat-1, -1), np.arange(min_lon, max_lon+1)
x_idxs, y_idxs = np.meshgrid(lon_range, lat_range)
y_meters = (y_idxs - min_lat) * dist_betn_lats
y_lats = y_idxs + min_lat
dist_betn_lons_at_min_lat = long_meters_at_lat(lat_range[-1])
x_meters = (x_idxs - central_lon) * dist_betn_lons_at_min_lat # Plus/minus around central longitude
x_lons = central_lon + np.round(x_meters/long_meters_at_lat(lat_range)[:, None]).astype('uint16')
assert ((x_lons[:, -1] - x_lons[:, 0]) <= 360).all(), 'The area is wrapping around on itself'
x_lons = np.where(x_lons >= 360, x_lons-360, x_lons)
这就是y_lats, x_lons的样子，看起来很正常（右上角的低经度绕了360°）。
(array([[75, 75, 75, ..., 75, 75, 75],
[74, 74, 74, ..., 74, 74, 74],
[73, 73, 73, ..., 73, 73, 73],
...,
[ 2, 2, 2, ..., 2, 2, 2],
[ 1, 1, 1, ..., 1, 1, 1],
[ 0, 0, 0, ..., 0, 0, 0]]),
array([[ 87, 91, 94, ..., 66, 69, 73],
[ 97, 101, 104, ..., 56, 59, 63],
[107, 110, 113, ..., 47, 50, 53],
...,
[215, 216, 217, ..., 303, 304, 305],
[215, 216, 217, ..., 303, 304, 305],
[215, 216, 217, ..., 303, 304, 305]], dtype=uint16))
2。我如何在地球上绘制这些 latitudes/longitudes 处的数据？
我在下面尝试了明显的，但只是在右边有一条窄条。
crs = ccrs.PlateCarree()
u = np.random.rand(*x_lons.shape)
v = np.random.rand(*x_lons.shape)
ax = plt.axes(projection=ccrs.Orthographic())
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax.set_global()
ax.scatter(y_lats, x_lons, u, v, transform=crs)
plt.show()
明显的错误是代码中 (long, lat) 的反转。 这是要尝试的正确代码。
# (second part only)
crs = ccrs.PlateCarree()
u = np.random.rand(*x_lons.shape)
v = np.random.rand(*x_lons.shape)
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-80, central_latitude=30))
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax.set_global()
ax.scatter(x_lons, y_lats, u, v, transform=crs)
plt.show()
编辑 1
这里是完整的代码，只在地图上的某个矩形形状内绘制数据。
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
import numpy as np
import matplotlib.patches as mpatches
# part 1 (minor change)
def long_meters_at_lat(lat):
"""Calculate distance (in meters) between longitudes at a given latitude."""
a = 6378137.0
b = 6356752.3142
e_sqr = a**2 / b**2 -1
lat = lat * 2 * np.pi / 360
return np.pi * a * np.cos(lat) / (180 * np.power(1 - e_sqr * np.square(np.sin(lat)), .5))
min_lat, max_lat = 0, 75
min_lon, max_lon = 215, 305 # Desired longitude spread at at min_lat
central_lon = (max_lon + min_lon) // 2
mean_lat = (max_lat + min_lat) // 2
dist_betn_lats = 111000 # In meters. Roughly constant
lat_range, lon_range = np.arange(max_lat, min_lat-1, -1), np.arange(min_lon, max_lon+1)
x_idxs, y_idxs = np.meshgrid(lon_range, lat_range)
y_meters = (y_idxs - min_lat) * dist_betn_lats
y_lats = y_idxs + min_lat
dist_betn_lons_at_min_lat = long_meters_at_lat(lat_range[-1])
x_meters = (x_idxs - central_lon) * dist_betn_lons_at_min_lat # Plus/minus around central longitude
x_lons = central_lon + np.round(x_meters/long_meters_at_lat(lat_range)[:, None]).astype('uint16')
assert ((x_lons[:, -1] - x_lons[:, 0]) <= 360).all(), 'The area is wrapping around on itself'
x_lons = np.where(x_lons >= 360, x_lons-360, x_lons)
# part 2
from_lonlat_degrees = ccrs.PlateCarree()
# map projection to use
proj1 = ccrs.Orthographic(central_longitude=central_lon, central_latitude=mean_lat)
u = np.random.rand(*x_lons.shape) # 0-1 values
v = np.random.rand(*x_lons.shape)
# auxillary axis for building a function (lonlat2gridxy)
axp = plt.axes( projection = proj1 )
axp.set_visible(False)
# this function does coord transformation
def lonlat2gridxy(axp, lon, lat):
return axp.projection.transform_point(lon, lat, ccrs.PlateCarree())
fig = plt.figure(figsize = (12, 16)) # set size as need
ax = plt.axes(projection=proj1)
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
# create rectangle for masking (adjust to one's need)
# here, lower-left corner is (-130, 15) in degrees
rex = mpatches.Rectangle( ax.projection.transform_point(-130, 15, ccrs.PlateCarree()), \
6500000, 4500000, \
facecolor="none")
ax.add_artist(rex)
bb = rex.get_bbox() # has .contains() for use later
# plot only lines (x,y), (u,v) if their
# (x,y) fall within the rectangle 'rex'
sc = 1. # scale for the vector sizes
for xi,yi,ui,vi in zip(x_lons, y_lats, u, v):
for xii,yii,uii,vii in zip(xi,yi,ui,vi):
xj, yj = lonlat2gridxy(axp, xii, yii)
# check only p1:(xj, yj), can also check p2:(xii+uii*sc, yii+vii*sc)
# if it is inside the rectangle, plot line(p1,p2) in red
if bb.contains(xj, yj):
ax.plot((xii, xii+uii*sc), \
(yii, yii+vii*sc), \
'r-,', \
transform=from_lonlat_degrees) #plot 2 point line
pass
# remove axp that occupies some figure area
axp.remove()
# without set_global, only rectangle part is plotted
ax.set_global() # plot full globe
plt.show()