Resample and plot satpy scene data from geos to mercator using pyresample

Question

I have a scene object from eumetsat satellite.

My area of interest: (xmin_corner, ymin_corner, xmax_corner, ymax_corner) [40,-10,105,40] degrees

looks like this

area definition looks like this:

Area ID: msg_seviri_unknown_3km
Description: MSG SEVIRI unknown area definition with 3 km resolution
Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '45.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1848
Number of rows: 1665
Area extent: (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914)

I want to project this into a mercator projection for flat viewing.

I create an area definition using the scene area details from above.

                my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
                            width=1848, height=1665,
                            area_extent= (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914))
                
                new_scn = scn.resample(my_area)
                
                new_scn.show('natural_color')

and check the new data.

projection looks fine, but some area is being clipped.

I so then try and specify area extent degrees explicitly

my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
            width=1848, height=1665,
            area_extent= [40,-10,105,40], units='degrees')

new_scn = scn.resample(my_area)

new_scn.show('natural_color')

Now I get a better image, but with some skew. its shrunk on the horizontal axis. I understand that this happens, because I have specified the width and the height?

How do I choose the best height, width for best optimal realistic output post projecting?

Now, I want to produce an image, maintaining the realistic looks (aspect ratio?,) when i move from geos to merc, without specifying some height and width and skewing the image. I want the image to look something like this

Additional information: I think because my dataset has np.nans, because its at the edge of the view of the satellite, this creates some issues.

When i convert my area extent from the original scene to lon lats.

scn_area = scn['natural_color'].attrs['area']
# xmin, ymin
scn_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (39.95708481758364, -10.007507617303142) - same as original. 40, -10

# xmax, ymax
scn_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (inf, inf)

Is this the intended behaviour when computing area extents? Should not i get a lon lat explicitly, from projection coordinates, and not inf?

I'm guessing this inf inf is because of the nan values at the corner, where in the geos projection its not earth, but space.

now. when I do the same with my new merc area definition I get this.

my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
            width=1848, height=1665,
            area_extent= (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914))

#xmin, ymin
my_area.get_lonlat_from_projection_coordinates(-604581.2379, -1096647.3571)
# (40.068954335025296, -9.867939236404165) - close to original (40, -10)

#xmax, ymax
my_area.get_lonlat_from_projection_coordinates(4940163.8125, 3899023.914)
# (89.87824658822916, 33.2040663665062) - completely off. and im guessing this is why the clipping happens.

I realise this must be something to do with how when the Scene object is initialised, the area extent is automatically calculated or something, because of the outer space np.nans. And by explicitly mentioning my area extent in degrees, I'm able to kind of get over this problem.

Now. Assuming what I've done and my understanding is correct. I want to be able to get the final image in its original aspect ratio (if i can call it that), or what is best practice used on google maps etc.

score 1 · Accepted Answer · answered Jun 21 '23 at 13:48

There are a lot of different variables at play here and you've covered quite a few scenarios in one single question. I'll try to answer what I can.

First, it isn't clear to me if you're trying to match some external geographic area definition with your area definition. So forgive me if there are variables here that you can't/don't want to change. In your first case with your extents defined in meters you said some things are "clipped". How so? If they are, then why not increase your extents? Your area has a resolution of ~3km x ~3km. You could remove the width/height from your call and instead specify resolution=(3000, 3000) and then modify your extents and let Pyresample figure out how many pixels there need to be.

In your second case with the degrees extents, this gets tricky. You can print the area and the pixel_size_x/y properties to see that the pixel resolution is much larger than your original 3km area and that they are not equal in the x and y dimensions:

In [6]: my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
   ...:             width=1848, height=1665,
   ...:             area_extent= [40,-10,105,40], units='degrees')

In [7]: my_area
Out[7]:
Area ID: my_area
Description: my_area
Projection: {'datum': 'WGS84', 'k': '1', 'lon_0': '45.5', 'no_defs': 'None', 'proj': 'merc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1848
Number of rows: 1665
Area extent: (-612257.1994, -1111475.1029, 6623509.7022, 4838471.3981)

In [8]: my_area.pixel_size_x
Out[8]: 3915.458280066441

In [9]: my_area.pixel_size_y
Out[9]: 3573.5414419900067

The odd arc of missing data on the right may be due to the "reduce_data" functionality in Satpy. This is an ongoing issue that we're trying to improve. You can work around it by passing reduce_data=False to your Scene.resample call and see if it improves that. If that doesn't fix it, let me know.

In your last set of questions, you say the mercator projection coordinates to lon/lat are "completely off". What do you mean by that? How so? You're right that the Infs in the geostationary area are from space pixels. In general, you have to be careful of what units you are using and in what coordinate system those coordinates are valid. Meters in the geostationary projection are not the same meters in the mercator projection. Converting to lon/lat degrees only confuses things further since the min/max x/y in the geostationary projection have nothing to do with the min/max lon/lat in degree space and likely won't be the min/max in mercator space. Since they are different coordinate systems you can't guarantee that the left-most part of the image in one is the left-most part of the image in another. Hopefully this rambling makes some sense.

My suggestion: If you know you want a 3km x 3km area in mercator and want it to cover a specific region in degree space then:

import xarray as xr
my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
    resolution=xr.DataArray([3000, 3000], attrs={"units": "meters"}),
    area_extent= [40,-10,105,40], units='degrees')

I think this should work. It basically says "extents are in degrees, resolution is in meters in the mercator projection".

My guess is reduce_data=False will also be needed to get around issues with cropping/clipping geostationary data.

Thanks for covering so much in one answer. setting the resolution using xarray and specifying units as meters did the trick for me. When i was trying to give resolution as 3000, earlier it was failing because units i was specifying degrees. — Rex, Jun 22 '23 at 05:50

Resample and plot satpy scene data from geos to mercator using pyresample

1 Answers1