Index Error when calling `boundary` with non full disk geos ara

[BENR0]

Index Error when calling `boundary` with non full disk geos ara

[djhoese]

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:

#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:

pyresample/pyresample/geometry.py

Lines 2785 to 2793 in 07382cc

	from shapely.geometry import Polygon
	geo_bbox = Polygon(np.vstack((x, y)).T)
	area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y)))
	intersection = area_bbox.intersection(geo_bbox)
	try:
	x, y = intersection.boundary.xy
	except NotImplementedError:
	return [], []
	return np.asanyarray(x[:-1]), np.asanyarray(y[:-1])

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).

[djhoese]

Oh or is this a shapely or numpy change that is causing this?

[BENR0]

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.

[djhoese]

This might have to be a @mraspaud answer (and @pnuu if he's available). Basically I see two possible solutions:

1. Do a 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.
2. Update the geostationary bounding box functions to not only use the number of points as a starting point for the full disk geostationary bounding polygon, but also make sure the subset results return that number of points.

[BENR0]

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.