地球儀上に「矩形」を形成するために緯度と経度を選択するにはどうすればよいですか？

Question

私は緯度経度グリッド上の値を持つデータセットを持っています。このデータセットから選択して、北アメリカのほぼ完全な「長方形」をプロットする必要があります。何かlike thisが、北米の上に置か：

1.どのように私は私の緯度と経度を選ぶのですか？

経度は極に向かって収束するので、北に向かってより長い経度が必要です。

ここに私のハッキーとおそらく不正確な試みがあります。私はこれがcartopyの1ライナーだと推測していますが、私はどのような変容を求めているのか分かりません。

私の四角形は、0度から75度までのN度の緯度です。私は、横幅（メートル）が緯度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) 

これは正気（右上の低経度は約360°をラップしている）思われるものをy_lats, x_lons表情のように、です。

(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.どのように私は地球上のこれらの緯度/経度のデータをプロットしますか？

私は以下で明らかにしようとしましたが、ちょうど右側に細いスライバがあります。

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() 

Answer 1

明白なエラーがあなたのコードで（、経度緯度）の逆です。 正しいコードを試してみてください。

# (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()