Open BENR0 opened 1 year ago
Oh no, what terrible timing to bring up this issue. @mraspaud just started his month long holiday and he has the most familiarity with this. It is possible some of this might get cleaned up with:
https://github.com/pytroll/pyresample/pull/526
since there was confusion caused by the naming of the frequency/number of pixels when generating the lon/lat arrays.
The frequency should have been set to 50
by default if I'm reading the source code correctly. So that's side02_step = 50 / 2 - 1 = 24
. So then the index here would be 24 * 2 + 1 = 49
, so that makes sense. The geostationary methods for generating the lons/lats seem to use np.linspace
to generate the original bounding arrays so then the only reason for this to get cut down to 7 (as the error message says) is if this intersection produces that small result:
So is it possible this intersection logic is not handling the anti-meridian or something? Oh but your example here is GOES East...that should be fine, right? You/we might need to add some prints to figure out why the geos boundary is so small (7 elements).
Oh or is this a shapely or numpy change that is causing this?
I saw #526 but did not read through it yet. From the code I was a little confused since frequency
did not really match up with the number of points I got which I guess is the reason for the mentioned PR. I don't think this is due to numpy or shapely since get_geostationary_bounding_box_in_lonlats
works fine. Additionally for a full disk AreaDefinition
there is no IndexError.
I think the reason is that for a non full disk AreaDefinition
the lat/lon arrays are calculated with the given frequency for the full disk but are later "cut" down due to the intersection but at the same time theboundary
method still uses the "full" frequency/number of points for the slicing.
This might have to be a @mraspaud answer (and @pnuu if he's available). Basically I see two possible solutions:
min
on the index used in the indexing operation shown in the original traceback. So it figures out the indexing it should used based on the frequency, but goes no further than the number of points in the array...is this possible? Maybe not because we don't know how to divide the polygon's coordinates into 4 sides which is what we're trying to do, right? Or...do we do min/max to get points that are on each half of the mid-lines (one horizontal midline, one vertical) to determine which points are left side, top side, etc.I think option 2. would be the cleaner one but might involve more refactoring. My first guess for just solving the problem was more on the line of option 1. if I understand it correctly.
When trying to get the boundary of a non full disk geostationary
AreaDefinition
currently anIndexError
will be thrown.Expected Output
The boundary lons/lats should be returned as in the full disk case.
Actual Result, Traceback if applicable
Versions of Python, package at hand and relevant dependencies
Python: 3.10.8 Pyresample: v1.27.1