djhoese
commented
7 months ago Try increasing your radius to 5000 and see how that changes the results.
Closed jpnIII closed 4 months ago
djhoese
commented
7 months ago Try increasing your radius to 5000 and see how that changes the results.
jpnIII
commented
7 months ago Thank you, Dave -- Please see my attached attempts with 5000 and 10000. No go. I am going to try to change my resolution/size parameters from: 1500 -1500 2000 1500 to: 1600 -1600 2000 1500
Thank you again, sir.
Jim
jpnIII
commented
7 months ago And here is my attempt with search radius = 5000, resolution/size parameters = 1600 -1600 2000 1500 I'm afraid this didn't work, either.
jpnIII
commented
7 months ago Thank you, Dave. I have tried your suggestion, as well as increasing to 10000. I also tried to change the resolution from 1500 to 1600. No go on all counts.
The weird thing is, everything worked beautifully with search radius = 1000, resolution/size = 1000 -1000 1600 1200
I think that the problems began when I went from size = 1600 1200 to 2000 1500 Same aspect ratio.
Would you suspect at all problems with odd integers? I would think not.
Later.
Sincerely,
Jim
From: David Hoese @.> Sent: Thursday, April 11, 2024 9:11 AM To: ssec/polar2grid @.> Cc: James P. Nelson @.>; Author @.> Subject: Re: [ssec/polar2grid] Problems running geo2grid.sh on RadM1 imagery from 08Apr2024? (Issue #691)
Try increasing your radius to 5000 and see how that changes the results.
— Reply to this email directly, view it on GitHubhttps://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2049784150__;Iw!!Mak6IKo!Je2_4CipPQ3gp8VppUspO3eXypfg4-xZ4BmImfS9nRfYttBjl00rZ7n7rvj4IOxnGdJjgLRhGcHYPSX_pTqwe2RF7Y3q-A$, or unsubscribehttps://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AARELF2SF2UNFA5CFN5DESDY42KXRAVCNFSM6AAAAABGBVLQB2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANBZG44DIMJVGA__;!!Mak6IKo!Je2_4CipPQ3gp8VppUspO3eXypfg4-xZ4BmImfS9nRfYttBjl00rZ7n7rvj4IOxnGdJjgLRhGcHYPSX_pTqwe2R_r9FNew$. You are receiving this because you authored the thread.Message ID: @.**@.>>
djhoese
commented
7 months ago So good news is that I can reproduce this with your data files. The bad news is I don't have much of an idea what is going on. At first I thought it was something with the LCC projection and maybe your grid was on the edge of it, but this doesn't make sense since your coastlines are shown fine. I then played with difference sizes and I think I've determined that it has something to do with the upper extent of your grid. I did a resolution of 1500x1500 and size of 1170Hx2000W and it was a bad image. If I change the center from your original 37 to 36, the image is fine.
I'll have to think about this, but at this point I'm not really sure why this upper limit is such a problem. It is possibly going beyond the Earth disk of the GOES-16 projection, but it doesn't make sense to me why that's a problem.
jpnIII
commented
7 months ago Thank you for checking into this, Dave. BTW, I tried checking into PROJ.4 projections and found a "US_Lambert_Conformal_Conic" projection. The navigation .yaml file generated was a bit different that the default, so I gave it a go within geo2grid.sh. Still did not work
.
FYI, attached is a GIF generated in McIDAS-X that shows the actual imagery coverage of the GOES-16 M1 sector image from 190027 UTC 08Apr2024 (099).
jpnIII
commented
7 months ago Thank you for checking into this, Dave. BTW, I tried checking into PROJ.4 projections and found a "US_Lambert_Conformal_Conic" projection. The navigation .yaml file generated was a bit different that the default, so I gave it a go within geo2grid.sh. However, this still did not work – similar to previous results.
I will give the Mercator projection a try.
Later.
Sincerely,
Jim
From: David Hoese @.> Sent: Thursday, April 11, 2024 3:49 PM To: ssec/polar2grid @.> Cc: James P. Nelson @.>; Author @.> Subject: Re: [ssec/polar2grid] Problems running geo2grid.sh on RadM1 imagery from 08Apr2024? (Issue #691)
So good news is that I can reproduce this with your data files. The bad news is I don't have much of an idea what is going on. At first I thought it was something with the LCC projection and maybe your grid was on the edge of it, but this doesn't make sense since your coastlines are shown fine. I then played with difference sizes and I think I've determined that it has something to do with the upper extent of your grid. I did a resolution of 1500x1500 and size of 1170Hx2000W and it was a bad image. If I change the center from your original 37 to 36, the image is fine.
I'll have to think about this, but at this point I'm not really sure why this upper limit is such a problem. It is possibly going beyond the Earth disk of the GOES-16 projection, but it doesn't make sense to me why that's a problem.
— Reply to this email directly, view it on GitHubhttps://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2050512619__;Iw!!Mak6IKo!JT_3oWMgZAALm3wW54-RqVHr_JylazEAGdg3bHGkXEhq3mxNJiVeUfvpvvR2iEeGcyquRZL9U5S1nNExlIR8ktL9i0A0iA$, or unsubscribehttps://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AARELFYQUIMUFQGHBVHI6OTY43ZKTAVCNFSM6AAAAABGBVLQB2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANJQGUYTENRRHE__;!!Mak6IKo!JT_3oWMgZAALm3wW54-RqVHr_JylazEAGdg3bHGkXEhq3mxNJiVeUfvpvvR2iEeGcyquRZL9U5S1nNExlIR8ktL0gfoPwQ$. You are receiving this because you authored the thread.Message ID: @.**@.>>
jpnIII
commented
7 months ago Hi again, Dave -- I thought I'd try a previous version of Geo2Grid, so I tried version 1.1 (installed on my computer back in Dec. 2022). Perhaps the following output will help you out:
. . . DEBUG : Warning: invalid value encountered in true_divide
DEBUG : Warning: invalid value encountered in true_divide
DEBUG : Warning: invalid value encountered in true_divide
WARNING : Resampling found -129405.02% of the output grid 'mesoscale_m1' covered. Will skip producing this product: C03
DEBUG : Unloading 'DataID(name='C03', wavelength=WavelengthRange(min=0.8455, central=0.865, max=0.8845, unit='µm'), resolution=1000, calibration=
djhoese
commented
7 months ago I've had time to look into this more and have narrowed it down. You/we are running into a long standing issue in satpy/pyresample related to a "reduce_data" option in resampling. Basically this functionality cuts off the parts of the input data that don't intersect with the output grid (AreaDefinition) to save processing time. Right now there is no way to disable this logic in Geo2Grid. For some reason, the intersection seems valid, but the resulting calculations are determining that the intersection only uses/needs 1 row of the input data. Obviously that isn't correct.
For my own record keeping as I debug this, the intersection is a box of:
# lons
[-99.32561454 -84.88516832 -83.0471571 -94.43288115]
# lats
[44.8036327 44.24231052 31.16295788 31.41858294]This maps to the following X/Y coordinates in the geostationary projection:
# x
[-1743495.04056 -741486.39656 -741486.39656 -1743495.04055999]
# y
[4188396.13192001 4188396.13192001 3186387.48792001 3186387.48791999]
# ABI area extent (minx, miny, maxx, maxy)
(-1743495.04056, 3186387.48792, -741486.39656, 4188396.13192)So the intersection of the input and output grid match the input data it seems, but perhaps some floating point precision is ruining it.
More debugging to come...
djhoese
commented
7 months ago @mraspaud @pnuu are the satpy and pyresample devs who have run into this issue in the past. I'm going to dump some findings here for them to provide feedback if they so desire:
The array indices computed from the above x/y coordinates end up being things like -0.50000001 and -0.499999998 or some floating point nonsense like that. There is also the other side of that for the 1000x1000 pixel area where the maximum index is 999.5000001 and 999.499999 or something close to that. The point is that the masked_ints decorator added to AreaDefinition.get_array_indices_from_lonlat which rounds the indices to integers and then masks anything < 0 and >= self.width/height. With floating point precision issues and with area extents being a half pixel beyond the array index, I don't think this decorator masked_ints works for this logic. I added a little blip to the code:
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 70e1c62..673886a 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -1447,10 +1447,13 @@ class _ProjectionDefinition(BaseDefinition):
def masked_ints(func):
"""Return masked integer arrays when returning array indices."""
@wraps(func)
- def wrapper(self, xm, ym):
+ def wrapper(self, xm, ym, allow_half_pixels=False):
is_scalar = np.isscalar(xm) and np.isscalar(ym)
x__, y__ = func(self, xm, ym)
+ if not is_scalar and allow_half_pixels:
+ x__ = np.clip(x__, 0, self.width - 1)
+ y__ = np.clip(y__, 0, self.height - 1)
x__ = np.round(x__).astype(int)
y__ = np.round(y__).astype(int)Then passed allow_half_pixels=True in pyresample/future/geometry.py line 57 in the get_area_slices method. I just did it this way because I didn't want to break anything. And in fact, this fixes my issue discovered here and gets me a full dataset. I'm not sure if this is something that needs to be considered for all masked_ints usage or just this area slices conditions. One more flexible option to my above hack would be to mask things that are outside the area array indices plus half a pixel plus some epsilon.
jpnIII
commented
7 months ago Thank you again for this work, Dave – You seem to be right on the cusp of solving the problem!
Sincerely,
Jim From: David Hoese @.> Sent: Thursday, April 18, 2024 8:31 PM To: ssec/polar2grid @.> Cc: James P. Nelson @.>; Author @.> Subject: Re: [ssec/polar2grid] Problems running geo2grid.sh on RadM1 imagery from 08Apr2024? (Issue #691)
@mraspaudhttps://urldefense.com/v3/__https:/github.com/mraspaud__;!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI6-ktrUDeQ$ @pnuuhttps://urldefense.com/v3/__https:/github.com/pnuu__;!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI682_WTmgg$ are the satpy and pyresample devs who have run into this issue in the past. I'm going to dump some findings here for them to provide feedback if they so desire:
The array indices computed from the above x/y coordinates end up being things like -0.50000001 and -0.499999998 or some floating point nonsense like that. There is also the other side of that for the 1000x1000 pixel area where the maximum index is 999.5000001 and 999.499999 or something close to that. The point is that the masked_ints decorator added to AreaDefinition.get_array_indices_from_lonlat which rounds the indices to integers and then masks anything < 0 and >= self.width/height. With floating point precision issues and with area extents being a half pixel beyond the array index, I don't think this decorator masked_ints works for this logic. I added a little blip to the code:
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 70e1c62..673886a 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -1447,10 +1447,13 @@ class _ProjectionDefinition(BaseDefinition):
def masked_ints(func):
"""Return masked integer arrays when returning array indices."""
@wraps(func)def wrapper(self, xm, ym):
def wrapper(self, xm, ym, allow_half_pixels=False):
is_scalar = np.isscalar(xm) and np.isscalar(ym)
x__, y__ = func(self, xm, ym)if not is_scalar and allow_half_pixels:
x = np.clip(x, 0, self.width - 1)
y = np.clip(y, 0, self.height - 1)
x__ = np.round(x__).astype(int)
y__ = np.round(y__).astype(int)Then passed allow_half_pixels=True in pyresample/future/geometry.py line 57 in the get_area_slices method. I just did it this way because I didn't want to break anything. And in fact, this fixes my issue discovered here and gets me a full dataset. I'm not sure if this is something that needs to be considered for all masked_ints usage or just this area slices conditions. One more flexible option to my above hack would be to mask things that are outside the area array indices plus half a pixel plus some epsilon.
— Reply to this email directly, view it on GitHubhttps://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2065585772__;Iw!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI6-17Qf4IA$, or unsubscribehttps://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/AARELF5LMA3GTZFWCEKLKP3Y6BXUJAVCNFSM6AAAAABGBVLQB2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANRVGU4DKNZXGI__;!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI6_ZXgMXJg$. You are receiving this because you authored the thread.Message ID: @.**@.>>
pnuu
commented
7 months ago diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 70e1c62..673886a 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -1447,10 +1447,13 @@ class _ProjectionDefinition(BaseDefinition): def masked_ints(func): """Return masked integer arrays when returning array indices.""" @wraps(func) - def wrapper(self, xm, ym): + def wrapper(self, xm, ym, allow_half_pixels=False): is_scalar = np.isscalar(xm) and np.isscalar(ym) x__, y__ = func(self, xm, ym) + if not is_scalar and allow_half_pixels: + x__ = np.clip(x__, 0, self.width - 1) + y__ = np.clip(y__, 0, self.height - 1) x__ = np.round(x__).astype(int) y__ = np.round(y__).astype(int)
My gut feeling is that the clipping should always be used. The decorator is used only twice and both seem to get array indices so where this issue will give unexpected results.
mraspaud
commented
7 months ago My understanding of this modification is that the clipping renders the masking useless, and thus no data is ever masked.
So if I provide a lon/lat array that produces eg [-1, 0, 1, 2, ...] for x and [..., 7, 8, 9, 10] for y (supposing a 10x10 array) then masked ints will not return [-, 0, 1, 2, ...] for x and [..., 7, 8, 9, -] for y, but [0, 0, 1, 2] for x and [7, 8, 9, 9] for y. I'm not sure what the effect would be if this is then used for indexing an array.
Or am I missing something?
djhoese
commented
7 months ago @mraspaud You are not missing anything, but @pnuu is right too that masked_ints is only used for the two "array indices". I think the main takeaway from exploring this is that there are two uses cases for these methods but only one is covered by the code currently:
get_area_slices logic. When we have a boundary instead of a 2D array of lon/lats (a "raster" of coordinates) then masking doesn't make sense. BUT maybe this clipping shouldn't be in these methods? I'm not sure.pnuu
commented
7 months ago Oh, right. I missed that the other use is for lonlats.
djhoese
commented
7 months ago Ok here's a slightly less controversial change I think:
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 70e1c62..0a0250b 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -1451,11 +1451,14 @@ def masked_ints(func):
is_scalar = np.isscalar(xm) and np.isscalar(ym)
x__, y__ = func(self, xm, ym)
+ epsilon = 0.02 # arbitrary buffer for floating point precision
+ x_mask = ((x__ < -0.5 - epsilon) | (x__ > self.width - 0.5 + epsilon))
+ y_mask = ((y__ < -0.5 - epsilon) | (y__ > self.height - 0.5 + epsilon))
+ x__ = np.clip(x__, 0, self.width - 1)
+ y__ = np.clip(y__, 0, self.height - 1)
x__ = np.round(x__).astype(int)
y__ = np.round(y__).astype(int)
- x_mask = ((x__ < 0) | (x__ >= self.width))
- y_mask = ((y__ < 0) | (y__ >= self.height))
x_masked = np.ma.masked_array(x__, mask=x_mask, copy=False)
y_masked = np.ma.masked_array(y__, mask=y_mask, copy=False)
if is_scalar:Plus sides: This fixes these index methods including up to 0.5 whole pixel more than they should on the right (the mask should have been self.width - 0.5 to include the outer extent but not further) and 0.5 less pixels on the left (the mask limit should have been -0.5). It includes some wiggle room by using 0.02 outside of the area extents to avoid floating point issues with lon/lat -> x/y and other calculations.
This has a downside that the masked pixels are now clipped/modified whereas before they remained their invalid out-of-bounds index values. I did it this way because it avoids me having to do weird complex clipping logic like x__[(x__ > -0.5 - epsilon) && (x__ < 0)] = 0.
mraspaud
commented
7 months ago Your last solution looks like a good compromise. However, maybe that function being used for two different purposes shows we should actually have two functions instead of one?
djhoese
commented
7 months ago That was going to be my other idea before I came up with the masked one above. However, I'm not sure these are technically different use cases, but rather stricter rules for how the methods/functions need to behave. I think there is a method that returns the floating point row/col indexes and if so then the get_area_slices could use that and do the rounding as needed. My above solution though definitely fixes some bugs though so I think it is still accurate. I could do more work to retain the invalid values underneath the mask if that is desired, but otherwise the above solution still fits in with the previous behavior.
djhoese
commented
3 months ago Oh interesting, I was able to close this issue from the PR in pyresample. FYI this is now released in pyresample 1.29.0, but I'll need to update the version used in the tarball as part of Geo2Grid.
Hi all -- I am having problems running geo2grid.sh on RadM1 sector GOES-16 ABI imagery from 08Apr2024 (099), trying to generate cimss_true_color and C03 imagery. I am just trying to get a C03 image now, since that is simpler than cimss_true_color. I am using Geo2Grid version 1.2.b. I am using ABI imagery from 19:00:27, and I know that all 16 bands exist for this date and time. I am using a center of (37.328,-89.777), which I also know is fine. My search radius is 1000. I am right in the middle of this mesoscale imagery, and should be able to use any search radius (?) and still get an image. Or so I think. The difficult thing is that I can get a nice cimss_true_color image using a resolution of (1000,-1000), size of (1600,1200), but then when I change to, for example, resolution = (1500,-1500), size = (2000,1500) I end up with a nonsense remapped image . The output says "SUCCESS", but the output image has problems. I have attached my output jpg (converted the tif from geo2grid.sh) file, as well as a png file with a map on it. Finally, I am trying to do this on imaginator.ssec.wisc.edu.
Thank you for your help!
Sincerely,
Jim
MY COMMANDLINE:
MY .yaml navigation file:
(make_geo2grid_files_with_paramfile16.scr) Following are all grid definitions within the file: geo2grid_config_mesoscale_m1.yaml:
MY LOGFILE:
MY OUTPUT FILE: