Resampling GOES mesoscale data to my area gives blank data

[agp-earth]

I am trying to resample GOES 17 mesoscale data. When I resample it to an area I defined with geometry.AreaDefinition, it crops but the data is all NaN. I tried the same code but with GOES 17 L1b Radiance data and that resamples to the area and outputs the correct data.

scn = Scene(filenames=glob.glob(myfile), reader='abi_l1b')
scn.load(['C07'],calibration='radiance')
projection = '+proj=utm +zone=3 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
proj_dict = utils.proj4_str_to_dict(projection)
area_def = geometry.AreaDefinition(area_id, description, proj_id, proj_dict, x_size, y_size, area_extent)
cropscn = scn.resample(my_area)
cropdata = cropscn['C07'].data
nan_check = np.isnan(np.array(cropdata.max()))

[djhoese]

Please provide the size and extent you used.

[agp-earth]

centerpt = np.array([6068126.51, 566285.63])
cellsize = scn['C07'].resolution
x_size = 100
y_size = 100
area_extent = ((centerpt[0]-(x_size/2.)*cellsize),(centerpt[1]-(y_size/2.)*cellsize),(centerpt[0]+(x_size/2.)*cellsize),(centerpt[1]+(y_size/2.)*cellsize))
proj_dict = utils.proj4_str_to_dict(projection)
area_def = geometry.AreaDefinition(area_id, description, proj_id, proj_dict, x_size, y_size, area_extent)

[agp-earth]

I am using this file: OR_ABI-L1b-RadM1-M6C07_G17_s20200080815309_e20200080815378_c20200080815417.nc which is a date where the mesoscale 1 overlaps with my area.

[djhoese]

Try passing radius_of_influence=50000 to the .resample call. My guess is this is so close to the edge of the geostationary disk that the default guessing isn't going to produce anything good.

Side note: There should be no need to use utils.proj4_str_to_dict anymore. Just pass the string directly to the AreaDefinition as the projection.

[agp-earth]

Thank you! I passed radius_of_influence=50000 into the resample call but it did not change anything. I also tried to make the area smaller. It does resample correctly if ,instead of using my area, I pass resampler='native', but then I cannot crop it to the desired area.

[djhoese]

I typed all of the below and then realized you said in your original post that the G17 radiance case was working for you? I'm confused now at what wasn't working. Anyway...

I tested a version of your original code with GOES-17 2022 data, but I used full disk because I had it on my machine already and I wanted to prove it wasn't some other part of the processing causing your issue.

So with full disk G17 data I see a proper image with:

centerpt = np.array([6068126.51, 566285.63])
cellsize = scn['C07'].resolution
x_size = 100
y_size = 100
area_extent = ((centerpt[0]-(x_size/2.)*cellsize),(centerpt[1]-(y_size/2.)*cellsize),(centerpt[0]+(x_size/2.)*cellsize),(centerpt[1]+(y_size/2.)*cellsize))
area_def = geometry.AreaDefinition("", "", "", projection, x_size, y_size, area_extent)

scn = Scene(reader="abi_l1b", filenames=glob("/data/satellite/abi/g17_2022061/OR_ABI-L1b-RadF-M6C07_G17_s20220611940321_e20220611949398_c20220611949428.nc"))
scn.load(["C07"])

cropscn = scn.resample(area_def)
cropscn["C07"].plot.imshow().figure.show()

I had looked up your file and I think you're right that that mesoscale on that date is in the location you specified (I think). My other guesses:

1. You're specifying calibration="radiance". Do you need that? If you don't, does it work without it?
2. I wonder if .resolution that you're using is somehow invalid for the version of Satpy you have or for the file you're using. I'm not sure that is possible given how the reader is configured but yeah...
3. The data might not actually be where we think it is.

[agp-earth]

The G17 radiance full disk works for me, it is just with the mesoscale that I run into this problem.

I removed the calibration="radiance" but that didn't change anything. I also checked the area extent to make sure I am within the image bounds.

If instead of using my defined area, I use resampler='native', it does work for Mesoscale.

[djhoese]

If you use pyproj and transform the input mesoscale extents to your target area definition's CRS you see the coordinates are way different than the extents you've specified:

In [19]: from pyproj import Transformer

In [20]: trans = Transformer.from_crs(src_area.crs, area_def.crs)

In [22]: trans.transform(src_area.area_extent[0], src_area.area_extent[1])
Out[22]: (705962.4980286012, 5109369.184016362)

In [23]: trans.transform(src_area.area_extent[2], src_area.area_extent[3])
Out[23]: (716941.0587718873, 7703772.140405564)

If I change the center point to something near the middle of these coordinates I get output when I do:

In [25]: centerpt = [710000.0, 6500000.0]
    ...: cellsize = scn['C07'].resolution
    ...: x_size = 100
    ...: y_size = 100
    ...: area_extent = ((centerpt[0]-(x_size/2.)*cellsize),(centerpt[1]-(y_size/2.)*cellsize),(centerpt[0]+(x_size/2.)*cellsize),(centerpt[1]+(y_size/2.)*cellsize))
    ...: area_def = geometry.AreaDefinition("", "", "", projection, x_size, y_size, area_extent)

In [26]: cropscn = scn.resample(area_def)

In [27]: cropscn["C07"].plot.imshow().figure.show()

[djhoese]

Let me know if this reveals something for you. I'm going to close this in the mean time.

[agp-earth]

Thank you so much!! This was incredibly helpful!