1

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)


ax.add_feature(cartopy.feature.LAND, edgecolor='black',alpha=0.5,zorder=1)

enter image description here

The four corners of the bounding box are not rectangular (which is fine), but for some reason the four corners of the scene are.

Method 2

Mercator plot

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

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

Now the scene is non-rectangular, but the corners of the bounding box are rectangular

enter image description here

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)

enter image description here

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!

0

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

enter image description here

1
  • 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 Mar 4 at 21:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.