Question about cropping resampled scenes

[lauracrews]

Thanks for this incredibly helpful set of tutorial notebooks and YouTube video of the SciPy workshop. After working through them, I still have a question about cropping/restricting an image to a given bounding box.

I’m working with Level 1b VIIRS imagery. My aim is to make true color images, resample them to the NSIDC Sea Ice Polar Stereographic North projection EPSG:3413 with resampled horizontal resolution of 750 m, and then crop all of the images to be constrained to the same geographic bounds prior to saving them as GeoTIFFs.

I am able to successfully create a true color Scene, resample it, and save it as a GeoTIFF. My problem is with cropping all of the images to have the same lat/lon bounds before saving the GeoTIFFs. Following the 04_resampling.ipynb notebook, I tried cropped_resampled_scene = resampled_scene.crop(ll_bbox=image_crop_limits) with image_crop_limits=(-160.0, 72.5, -115, 80). However, I get the error “NotImplementedError: Source projection must be 'geos' if source/target projections are not equal” which makes me think crop only works for geostationary imagery prior to resampling?

In other words, I am currently able to generate images that look like this:

But I would like images that look like this:

It is important that I retain the projection and horizontal resolution through the cropping.

The code I am using to create an example uncropped image is:

#Imagery file names
data_path = '~/test_images/'
filenames = glob(data_path+'*nc')
filenames

#Moderate and imagery resolution files and geolocation files for the example image
#I don’t seem to be able to upload these to github
#[VJ102IMG.A2019207.2254.021.2021053085916.nc,
VJ103IMG.A2019207.2254.021.2021053075151.nc,
VJ102MOD.A2019207.2254.021.2021053085916.nc,
VJ103MOD.A2019207.2254.021.2021053075151.nc]

#Setup the projection
area_id = 'nsidc_polar_north'
description = 'NSIDC Sea Ice Polar Stereographic North'
projection = 'EPSG:3413'
resolution = 750 #in meters
my_dynamic_area = create_area_def(area_id=area_id, description=description, projection=projection, resolution=resolution)

#Make the true color scene
scn = Scene(filenames=filenames, reader='viirs_l1b')
channel = 'true_color'
scn.load([channel], generate=False)
resampled_scene = scn.resample(my_dynamic_area)

#Crop the scene. This does not work
image_crop_limits=(-160.0, 72.5, -115, 80) #geographic coordinates to bound the final image
cropped_resampled_scene = resampled_scene.crop(ll_bbox=image_crop_limits)

#Save the image
resampled_scene.save_datasets()

Thanks so much for your help!
Laura

[djhoese]

Thanks for filing the issue. Could you provide the full traceback for that NotImplementedError you are seeing? I'm having trouble tracking down where in the code that is occurring. If you are using old versions of satpy and pyresample (like the ones originally used in this tutorial) then I highly recommend upgrading.

Otherwise, for your particular use case, it would be better to resample to a fully static (not dynamic) area with the exact crop limits you want. That is, unless you need to do some work with the "full swath" version of the resampled data. You could define your area definition with extents in lon/lat and let pyresample/satpy transform them to the proper values for your projection. Right now it requires some extra work to specify it this way but it should look like this:

import xarray as xr

ll_extents = xr.DataArray([-160.0, 72.5, -115.0, 80.0], attrs={"units": "degrees"})
my_area = create_area_def(area_id=area_id, description=description, projection=projection, resolution=resolution, area_extent=ll_extents)
# pass to .resample and don't .crop anymore

Reasons why this is better:

1. You remove the need for a crop call.
2. You avoid computing/freezing the coordinates for the dynamic area you were creating.
3. You end up resampling/producing less data.

[lauracrews]

Hi David,

Thanks so much for getting back to me, I really appreciate it. I am currently using Satpy 0.43.0. I see there was a recent release - happy to update if you think that's contributing.

Here's the full traceback for the error I was having with the dynamic area I was using:

However, it does sound like the static area with specified crop limits would be better for my case. I don't plan to do anything with the full swath data - rather, I plan to look at movement of ice floes over a series of images within the specified geographic region. Because I'll be tracking features over multiple images, I need to retain the same bounds and horizontal resolution.

When I tried implementing the proposed solution, I ran into new errors. Here is the traceback for the current warnings and errors:

Thanks again for your help with this! I've really appreciated how well supported and documented SatPy is.

[djhoese]

Ah it looks like the ability to do this was added right after the latest pyresample release and we haven't had one since:

pytroll/pyresample#516

Anyway, for your new error...could you print out my_area and show me the output. I have a feeling the create_area_def function is getting confused by the polar-centric projection you're using. If so, we can do this manually...I'll try doing the calculations myself while I wait for a reply.

[lauracrews]

Yes, absolutely. Here's the output of my_area:

[djhoese]

Yeah a negative number of rows is not a good sign.

So I did some of the projection transformation manually:

In [15]: from pyproj import Transformer

In [16]: transformer = Transformer.from_crs(4326, 3413)

In [17]: transformer.transform(-160, 72.5)
Out[17]: (inf, inf)

In [18]: transformer.transform(-160, 80)
Out[18]: (inf, inf)

In [19]: transformer = Transformer.from_crs(4326, 3413, always_xy=True)

In [20]: transformer.transform(-160, 72.5)
Out[20]: (-1731055.8346230083, 807204.5923786014)

In [21]: transformer.transform(-160, 80)
Out[21]: (-984178.0216285206, 458929.7484732148)

In [22]: transformer.transform(-115.0, 72.5)
Out[22]: (-1794821.1603503, -653261.4781986361)

In [23]: transformer.transform(-115.0, 80)
Out[23]: (-1020431.2902219342, -371406.61575464066)

We can see from the last set of calls that the minimum X coordinate is -1794821.1603503, the max X is -984178.0216285206, min Y is -653261.4781986361, and max Y is 807204.5923786014. If you plug those into extents you can do:

extents = (-1794821.1603503, -653261.4781986361, -984178.0216285206, 807204.5923786014)
my_area = create_area_def(..., area_extent=extents)

Then hopefully this should work. My guess is that the extents you chose, given the odd transform between lon/lat spaces and polar-stereographic, are not actually the "lower-left" and "upper-right" corner that pyresample assumed they were and things got weird after that.

[lauracrews]

That does seem to work! Thank you so much!

It would be really good if I understood your fix well enough to use it in the future, as it's likely I'll need other geographic ranges at some point. Does it generalize to something like:

from pyproj import Transformer
minlon = -160
maxlon = -115
minlat = 72.5
maxlat = 80

transformer = Transformer.from_crs(4326, 3413, always_xy=True)

#Define 'corners'
(x1, y1) = transformer.transform(minlon, minlat)
(x2, y2) = transformer.transform(minlon, maxlat)
(x3, y3) = transformer.transform(maxlon, minlat)
(x4, y4) = transformer.transform(maxlon, maxlat)

x = (x1, x2, x3, x4)
y = (y1, y2, y3, y4)

extent = (min(x), min(y), max(x), max(y))

[lauracrews]

One additional question: Is it possible to save the images without the alpha channel, as RGB instead of RGBA GeoTIFFs? I see there is the function exclude_alpha but I am not sure how to apply it. I tried the following:

Thank you!

[djhoese]

Does it generalize to something like:

Yes. Normally though this type of work isn't needed, but your polar projection/area is one of the more odd cases that pyresample tries to handle. I'm not familiar with this EPSG code or with the locations you specified but theoretically if you would have picked the lon/lats that corresponded to the lower-left and upper-right corners of the area in my original suggestion (with the DataArray and units metadata) it should have worked...but maybe that is what you did and it still didn't work. I'm not entirely sure. It could be a bug. It is at least some what of a small bug as it should have been more graceful with it's bad result.

Regarding no alpha: Pass fill_value=0 to save_datasets. The alpha channel in satpy is used to hide invalid values (transparent pixels for fill/NaN values) so for non-alpha versions of the image we need to Satpy (trollimage, the package doing the saving) what to do instead of this alpha channel generation. Passing fill_value=0 means "set invalid values to black".

[lauracrews]

Thanks for the explanations on both points! That makes sense. I really appreciate your assistance with these questions.