Plotting high latitude bounding box using cartopy

I have used the sentinelhub python api to download a Sentinel-1 scene corresponding to the north coast of greenland where a bounding box defined by lon-lat is highly distorted. I have tried a few things to plot it but none seem to match the cartopy land mask that I'm using.

``````# Code to get lon/lat coords

# This is adapted from this guide
# https://forum.sentinel-hub.com/t/get-output-coordinates-using-sentinelhubpy/2804/4

bbox = BBox([-40, 80, -5, 84], crs=CRS.WGS84)
bbox_size = bbox_to_dimensions(bbox, resolution=500)

bb_utm = sentinelhub.geo_utils.to_utm_bbox(bbox)
transf = bb_utm.get_transform_vector(resx=500, resy=500)

pix_row = np.arange(0, bbox_size[1])
rows = np.array([pix_row] * bbox_size[0]).transpose()

pix_col = np.arange(0, bbox_size[0])
cols = np.array([pix_col] * bbox_size[1])

# Convert the pixel positions to UTM
utmx, utmy = geo_utils.pixel_to_utm(rows, cols, transf)

# Convert your UTM pixel positions to WGS84 (EPSG:4326)
lon_degrees, lat_degrees = geo_utils.to_wgs84(utmx, utmy, bb_utm.crs)
``````

I have then tried several different ways of plotting this.

Method 1

North Polar Stereo with lon/lat coords

``````fig = plt.figure(figsize=(10,10))

ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=-30,true_scale_latitude=82))

# This is my bounding box
# [-40, 80, -5, 84]

ax.set_extent([-45, 0, 89, 79], ccrs.PlateCarree())

ax.pcolormesh(lon_degrees,
lat_degrees,
dtp[:-1,:-1],
vmin = 0,
vmax = 1,
transform=ccrs.PlateCarree(),
cmap='PRGn',
alpha=1)

ax.scatter([-40,-5],[80,84],transform=ccrs.PlateCarree(),color='r',zorder=2)
ax.scatter([-40,-5],[84,80],transform=ccrs.PlateCarree(),color='r',zorder=2)

``````

Method 2

Mercator plot

Here I've replaced the `plt.axes` call with:

``````ax = plt.axes(projection=ccrs.Mercator())
``````

Method 3

Do it all in UTM

``````ax = plt.axes(projection=ccrs.UTM(28))

# This is my bounding box
# [-40, 80, -5, 84]

ax.set_extent([-45, 0, 89, 79], ccrs.PlateCarree())

ax.pcolormesh(utmx,
utmy,
dtp[:-1,:-1],
vmin = 0,
vmax = 1,
transform=ccrs.UTM(28),
cmap='PRGn',
alpha=1)

``````

Bounding box coordinates are now non-rectangular again, and scene is rectangular.

Bizarrely this is the closest fit between image and land mask, but the top right and bottom left corners of the bounding box no longer correspond to the corners of the scene!

Okay so it turns out there is a considerably better match if you define your bounding box in UTM coordinates (although this is a bodge as you have to pick one zone that only covers a small area at high latitudes).

There are still some residual errors I think (visible near the top of the image), but it's a start.

I've put the code to produce the figure below here

• For answers that involve code we ask that you include that code (or at least a code snippet based on it) in your answer. There is an edit button beneath your question which will enable you to do that and a `{}` button that enables you to format any highlighted code nicely.
– PolyGeo
Commented Mar 4, 2021 at 21:44