Polar2Grid 3.1 not calculating grid coverage correctly

[kathys]

Douglas Schumacher submitted this P2G bug report to us:

I’ve been working with Liam to update packages on our processing system at HCC. When Polar2Grid was upgraded from v3.0 to v3.1, a problem started where, when CSPP MiRS v4.0 input is used, the grid-coverage attribute does not appear to be respected, resulting in GeoTIFFs being wrongly created where no data coverage exists. This behavior is consistent among all CSPP MiRS v4.0 input, but does not occur in any other input that I’ve seen (MODIS, VIRRS, AVHRR Level 1, for example). From what I can tell, P2G v3.0 handles the MiRS v4.0 files as expected (other than no NOAA-21 support), not producing output when there is no data coverage on the specified grid.

Here’s the command being used:
polar2grid.sh -r mirs -w geotiff -p btemp_31v btemp_88v btemp_165h btemp_183h1 btemp_183h5 rain_rate tpw --fill-value 0 -g --grid-configs storm_grid.conf --grid-coverage=0.45 --num-workers 8 -f /IMGnpp*.nc /IMGn2*.nc

Interestingly, in the logs, you can see P2G v3.1 recognizes that data from this pass doesn’t cover enough of the Honolulu grid to produce output, but for the other grids with zero data coverage, it decides to produce output:

INFO : Sorting and reading input files...
WARNING : Don't know how to open the following files: {'mirs/IMGnpp*.nc'}
INFO : Loading product metadata from files...
INFO : Loading swath geolocation into memory...
INFO : Running day coverage filtering...
INFO : Running night coverage filtering...
INFO : Checking products for sufficient output grid coverage (grid: '99B_INVEST')...
INFO : Resampling to '99B_INVEST' using 'ewa' resampling...
INFO : Checking products for sufficient output grid coverage (grid: '02S_INVEST')...
INFO : Resampling to '02S_INVEST' using 'ewa' resampling...
INFO : Checking products for sufficient output grid coverage (grid: '96S_INVEST')...
INFO : Resampling to '96S_INVEST' using 'ewa' resampling...
INFO : Checking products for sufficient output grid coverage (grid: '20E_INVEST')...
INFO : Resampling to '20E_INVEST' using 'ewa' resampling...
INFO : Checking products for sufficient output grid coverage (grid: 'HONOLULU')...
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: RR
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: TPW
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: btemp_165h
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: btemp_183h1
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: btemp_183h5
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: btemp_31v
WARNING : Resampling found 39.59% of the output grid 'HONOLULU' covered. Will skip producing this product: btemp_88v
WARNING : No products were found to overlap with 'HONOLULU' grid.
INFO : Computing products and saving data to writers...
INFO : SUCCESS

I’ve attached a sample of what one of these empty GeoTIFFs looks like – the actual grid is correct (coastlines and lat/lon in the right place), it just shouldn’t been generated because it doesn’t overlap with the MiRS input. Please let me know if you need any additional information to look into this. Thank you.

[kathys]

I requested this MiRS dataset from Douglas. I received it, and was able to reproduce the problem. The log files for the above image, which is the grid defined as 99B_INVEST, does not intersect with the actual data, yet the log file shows 100% coverage. I checked the MiRS input file and noted two things. The data is clean, with no missing data that I could see, and it cross the anti-meridian. Douglas also states that he does not remember this happening with P2G v3.0, but I could not test on this dataset because v3.0 did not support NOAA21.

I put the data and execution results here: bumi:/data/users/kathys/test_data/mirs/douglas_2024

[djhoese]

Just a small update on this: This is definitely the navigation of the swath crossing the anti-meridian that's the problem. Well, pyresample's inability to handle it properly I should say. It doesn't fall into the situation that we've seen with other data cases where the area of the bounding polygon is a negative number; we're still checking for that. If I add some code to Polar2Grid to ignore the bounding lon/lats of the swath and instead do a min/max of the entire swath then I get the correct result for 99B_INVEST of 0% coverage.

I/we know there are issues with pyresample on handling antimeridian and in the polygon intersection code in general so I'm not sure diving into the specifics of the current pyresample implementation is worth my/our time. I will think about if there is a short term check I can do to determine if I should do the min/max version of the bounding box instead of the "smarter" bounding lon/lats case.

[djhoese]

@kathys So I've come up with something, but it means slower processing speeds for every swath when grid coverage needs to be calculated. The solution is: if the longitudes of the bounding polygon have values above 170 and below -170, then we should consider all of the coordinates and make a min/max bounding box. So for small mirs swath where we're only looking at the outer edge pixels' coordinates this is nothing. For long I-band VIIRS swaths though this could be a couple seconds of extra time. I can't imagine it would be more than 10 seconds, but this is just me guessing.

[kathys]

I did not get any reports of VIIRS or any other products having the same issues. Douglas said it was only MiRS that this seemed to affect. Is that possible?

[djhoese]

Good point, but I think that might be more luck than something instrument-specific. It could be the spacing between the ATMS pixels that confuse the polygon math just enough to get this result.

[kathys]

I believe in luck! If we do not have to slow down the reprojections for other instruments/products, I would prefer that solution.

[kathys]

I should have added in the short term. Seems like there are some bigger things in pyresample that eventually need to be worked out.

[djhoese]

So it's pretty difficult to do this as instrument-specific due to where the problem code is. It only receives the lon/lats and makes a polygon for it. I did another quick test after remembering that the polygons are actually reduced from the full list of lon/lat side coordinates. I thought maybe if I didn't do the reduction maybe MiRS would do better, but no it still fails. I also tried reducing it more and that also didn't fix it. I'm not sure right now what is so terrible about these coordinates that the math doesn't like it.

[kathys]

I think your solution sounds like a good one for the time being Dave. Even using data from the Honolulu antenna, this is not going to happen for every overpass, and for the majority of other users, it will never be a problem.

[djhoese]

Plot of lon v lat for this case.

If I move all negative lons to the positive side (+ 360), much less sinister:

[djhoese]

I haven't come up with a better solution, but I did try more things. Even when I limit the polygon to the corners of the swath it still fails (100% coverage). The bbox solution which also has 4 corners passes, but the difference is it has two corners in the negative side and two on the positive side of longitudes. Something about the math must not like how the polygon switches for just one corner.