How do I select latitudes and longitudes to form a "rectangle" on a globe?

Question

I have a data set with values on a latitude-longitude grid. I need to select from this data set to plot a near-perfect "rectangle" over North America. Something like this, but placed over North America:

1. How do I choose my latitudes and longitudes?

Since longitudes converge towards the poles, I need more longitudes towards the north.

Here is my hacky and probably incorrect attempt. I'm guessing this is a one-liner in cartopy, but I don't know what transformation I'm looking for.

My rectangle has a height from 0° to 75°N latitudes. I'm calculating the span of longitudes at each latitude such that the horizontal width (in meters) is the same as the distance from 215° to 305° longitude at 0° latitude. (The rectangle is centered horizontally around 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)

This is what y_lats, x_lons look like, which seems sane (the low longitudes on the top right have wrapped around 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. How would I plot the data at these latitudes/longitudes on a globe?

I tried obvious below, but just get a narrow sliver off to the right.

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

The obvious error is the reverse of (long, lat) in your code. Here is the correct code to try.

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

Edit 1

Here is the complete code that only plot the data within a certain rectangular shape on the map.

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