Weird problem in `pyresample` newer version for the `VIIRS` data

[anikfal]

I'm using pyresample by the same code and on the same files, but in two separate workflows:

Code:

from satpy.scene import Scene

main375 = "VJ102IMG_NRT.A2023327.0830.021.nc"
geo375  = "VJ103IMG_NRT.A2023327.0830.021.nc"
main750 = "VJ102MOD_NRT.A2023327.0830.021.nc"
geo750  = "VJ103MOD_NRT.A2023327.0830.021.nc"
scn = Scene(reader='viirs_l1b', filenames=[main375, main750, geo375, geo750])
scn.load(["I01"])
my_area = scn["I01"].attrs['area'].compute_optimal_bb_area({'proj': 'lcc', 'lon_0': 35., 'lat_0': 50., 'lat_1': 50., 'lat_2': 50.})

Workflow 1 (no problem without errors): OS: CentOS 7 satpy version: 0.28.0 pyresample version: 1.19.0

Workflow 2 (with errors): OS: CentOS 7 satpy version: 0.43.0.post0 pyresample version: 1.27.1

Errors for workflow 2:

/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/xarray/core/dataset.py:270: UserWarning: The specified chunks separate the stored chunks along dimension "number_of_pixels" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/xarray/core/dataset.py:270: UserWarning: The specified chunks separate the stored chunks along dimension "number_of_pixels" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/xarray/core/dataset.py:270: UserWarning: The specified chunks separate the stored chunks along dimension "number_of_pixels" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/xarray/core/dataset.py:270: UserWarning: The specified chunks separate the stored chunks along dimension "number_of_LUT_values" starting at index 4096. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
Traceback (most recent call last):
  File "/home/sebal/app/viirs_make/zzz/reserve/band_maker.py", line 9, in <module>
    my_area = scn["I01"].attrs['area'].compute_optimal_bb_area({'proj': 'lcc', 'lon_0': 35., 'lat_0': 50., 'lat_1': 50., 'lat_2': 50.})
  File "/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/pyresample/geometry.py", line 990, in compute_optimal_bb_area
    lons, lats = self.get_edge_lonlats()
  File "/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/pyresample/geometry.py", line 420, in get_edge_lonlats
    lons, lats = self.get_bbox_lonlats(frequency=frequency, force_clockwise=False)
  File "/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/pyresample/geometry.py", line 331, in get_bbox_lonlats
    lons, lats = self._get_bbox_elements(self.get_lonlats, frequency)
  File "/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/pyresample/geometry.py", line 350, in _get_bbox_elements
    clean_dim1, clean_dim2 = self._filter_bbox_nans(dim1, dim2)
  File "/home/sebal/miniconda3/envs/mysat/lib/python3.10/site-packages/pyresample/geometry.py", line 363, in _filter_bbox_nans
    raise ValueError("Can't compute swath bounding coordinates. At least one side is completely invalid.")
ValueError: Can't compute swath bounding coordinates. At least one side is completely invalid.

---

This problem occurs occasionally. Perhaps once per month or even less. if I apply the code on some other VIIRS data, there won't be any errors. For example the files below:

VJ102IMG_NRT.A2023327.1012.021.nc
VJ102MOD_NRT.A2023327.1012.021.nc
VJ103IMG_NRT.A2023327.1012.021.nc
VJ103MOD_NRT.A2023327.1012.021.nc

[djhoese]

Thanks for filing the bug. Obviously a lot has changed between the versions of the libraries you're using (thanks for including that version information in your description).

In the old workflow, could you print out what scn["I01"].attrs["area"].get_edge_lonlats() and paste it here? My guess (hope?) is that although it maybe returned something that was used by the compute_optimal_bb_area that it was actually producing an inefficient/incorrect/inaccurate result. That would at least justify why the new version errors out and doesn't produce anything: it knows that the result isn't good or accurate.

NaNs in VIIRS or MODIS or other polar-orbiters is very common, especially on the first or last granule of a satellite pass. Unfortunately there usually/currently isn't a good way to deal with them while also having decent performance (not computing/loading every lon/lat and doing a min/max on them). I'll take a look at the methods involved after I have more information from you.

[anikfal]

You're right. Thanks. The data seems to be corrupted: scn["I01"].attrs["area"].get_edge_lonlats()

(masked_array(data=[nan, nan, nan, ..., nan, nan, nan], mask=False, fill_value=1e+20, dtype=float32), masked_array(data=[nan, nan, nan, ..., nan, nan, nan], mask=False, fill_value=1e+20, dtype=float32))

This is the first time I come across to a corrupted VIIRS data. However, that was interesting that pyresample v1.19.0 could handle or skip it, so the workflow was not terminated.

[djhoese]

Could you print out my_area for me? That is, the result of the omerc operations? Does the final image for that operation look OK?

[anikfal]

Sure. I insert the output image displayed by QGIS. pyresample v1.19.0 could at least make a GeoTIFF for the valid pixels (upper rows). The image, if was not corrupted, was supposed to cover eastern regions of Iran.

my_area = scn["I01"].attrs['area'].compute_optimal_bb_area({'proj': 'lcc', 'lon_0': 35., 'lat_0': 50., 'lat_1': 50., 'lat_2': 50.})

/home/sebal/miniconda3/envs/mysat/lib/python3.9/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj_string = self.to_proj4()

new_scn = scn.resample(my_area)
new_scn.save_datasets(writer='geotiff', dtype=np.float32, enhance=False)

[djhoese]

Sorry, I meant print(my_area).

[anikfal]

print(my_area)

/home/sebal/miniconda3/envs/mysat/lib/python3.9/site-packages/pyproj/crs/crs.py:543: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj_string = self.to_proj4()
Area ID: lcc_otf
Description: On-the-fly lcc area
Projection: {'ellps': 'WGS84', 'lat_0': '50', 'lat_1': '50', 'lat_2': '50', 'lon_0': '35', 'no_defs': 'None', 'proj': 'lcc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 16936
Number of rows: 6969
Area extent: (89458.8947, -725329.4873, 3348658.0572, 1987300.0824)

[djhoese]

Ok I had to skim through some of this code because I was starting to get confused why you were calling the compute_optimal_bb_area method directly since I'm used to calling .freeze, but anyway...

If we go back to lons, lats = scn["I01"].attrs["area"].get_edge_lonlats(), could you take that result and do np.nanmin and np.nanmax on lons.data and lats.data? Basically I'm curious if there are any non-NaNs in the result and what our limits are.

I think I still stand by the newer versions "opinion"/decision that producing a result from this many NaNs isn't going to be accurate. That said, I have implemented a workaround to basically this exact situation in my own project that uses pyresample that may I should implement here if it isn't already. Basically if we can't compute/freeze the dynamic area with the lon/lat boundaries, then do the costly computation of all the lon/lats and do a min/max on them. I'll try to look at this more this week.

[anikfal]

Here are the extra information for the corrupted VIIRS data, by pyresample v1.19.0. Let me know if more details are required.

lons, lats = scn["I01"].attrs["area"].get_edge_lonlats()
print(lons), print(lats)

[nan nan nan ... nan nan nan] [nan nan nan ... nan nan nan]

print( lons.shape ), print( lats.shape )

(25472,) (25472,)

np.all( np.isnan(lons) ), np.all( np.isnan(lats) )

(False, False)

sum( np.isnan(lons) ), sum( np.isnan(lats) ) #12352 nan-pixels out of 25472 pixels

(12352, 12352)

np.nanmax(lons), np.nanmin(lons)

(84.68566, 36.313953)

---

my_area = scn["I01"].attrs['area'].compute_optimal_bb_area({'proj': 'lcc', 'lon_0': 35., 'lat_0': 50., 'lat_1': 50., 'lat_2': 50.})
lons_optimal, lats_optimal = my_area.get_lonlats()
print( lons_optimal.shape ), print( lats_optimal.shape )

(6969, 16936) (6969, 16936)

np.any( np.isnan(lons_optimal) ), np.any( np.isnan(lats_optimal) ) #No nan-value for any array element

(False, False)

print( lons_optimal )

[[36.98399656 36.9882579 36.99251922 ... 93.44126797 93.4434173 93.44556652] [36.9837679 36.98802875 36.99228958 ... 93.43695573 93.43910506 93.44125427] [36.98353929 36.98779965 36.99205999 ... 93.43264399 93.43479332 93.43694253] ... [36.10037132 36.10273553 36.10509974 ... 72.61398324 72.61579863 72.61761398] [36.10030097 36.10266502 36.10502908 ... 72.61196351 72.61377884 72.61559412] [36.10023062 36.10259453 36.10495844 ... 72.60994397 72.61175924 72.61357446]]

np.max(lons_optimal) , np.min(lons_optimal)

(93.44556651556852, 36.10023062183153)

np.max(lats_optimal) , np.min(lats_optimal)

(67.51536906813563, 35.865271984165695)

[djhoese]

Ah very nice doing the optimal. We can see that the edge lonlats aren't including quite as much as the full array of lonlats would have. I don't see the edge lats min/max in your output but that's probably OK.

I did want to go over why this error is happening and what your code is doing even though I'm sure you know a lot of this information. This is partly for future me:

So the method you're calling is only being given the swath def (self in the internals of the method) and the projection you've provided. It is trying to determine what the best extents and number of rows and columns are for an area definition for that provided projection. So it first computes a reasonable number of rows and columns and since you aren't providing a preferred pixel resolution (this may not have been available in your old version of pyresample) it is trying to keep about the same number of rows and columns as the swath data. It also tries to add any projection parameters that you don't specify like lon_0 and lat_0 in your projection definition so they are related to the lon/lats of your swath definition. But you provide both of these parameters so the only thing I think it adds is the +ellps=WGS84 which is likely the default anyway. Lastly, it calls DynamicAreaDefinition.freeze with the computed rows/cols and the edge lon/lats which means it is essentially coming up with the "extents" of the area (or you could see it as determining the resolution). It is using these edge lon/lats versus the whole lon/lat array to save on performance. However, the edge lon/lats method doesn't know this so it is error'ing out because it doesn't make sense to have "edge" or boundary information if an entire side of that information is invalid. So while I think we could provide a workaround for the compute_optimal_bb_area I don't think the edge lon/lats method is where this should happen. CC @mraspaud who wrote most if not all of this code for optimal areas and should weigh in on this.

All of that said, I don't think you should be using this compute_optimal_bb_area method at all. For one, it gives you non-square pixels:

In [1]: extents = (89458.8947, -725329.4873, 3348658.0572, 1987300.0824)

In [2]: rows = 6969

In [3]: cols = 16936

In [4]: (extents[2] - extents[0]) / cols
Out[4]: 192.44208564596127

In [5]: (extents[3] - extents[1]) / rows
Out[5]: 389.2422972736404

Second, it isn't updating your projection at all as you already specify lon_0 and lat_0.

Third, the compute_optimal_bb_area method accesses lon/lats multiple times to compute those projection parameters and the rows and columns and then also the edge lon/lats. This is very wasteful and isn't helping performance much.

If you instead could determine a pixel resolution that would work for your data (375m for VIIRS I-band resolution?) and created a DynamicAreaDefinition with your projection and this resolution, then passed that to scn.resample you'd get essentially the same result with square pixels (same resolution in X and Y dimensions) and more consistent results...I think. Note that scn.resample is calling DynamicAreaDefinition.freeze for you. Again, @mraspaud may have different opinions on this and I may be missing something.

[mraspaud]

I agree with your analysis @djhoese and would also recommend not using compute_optimal_bb_area directly. For example, you could use an area like this:

laea_bb:
  description: Lambert Azimuthal Equal-Area Bounding Box for Polar Overpasses
  projection:
    ellps: WGS84
    proj: laea
  resolution: 300
  optimize_projection: True

(note the 300 resolution)

However, whether this will solve your problem or not, I don't know. Having NaNs at the edges is indeed a problem that we need to tackle, but as @djhoese, it's problematic performance-wise. A long term solution I believe would be to use provided metadata from the file (in viirs SDR files for example they to have a "geo ring" attribute that just give the granule edge lon-lats), but that's not on my short-term todo list yet.