Problems running geo2grid.sh on RadM1 imagery from 08Apr2024? #691
Comments
|
Try increasing your radius to 5000 and see how that changes the results. |
|
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
|
|
And here is my attempt with search radius = 5000, resolution/size parameters = 1600 -1600 2000 1500 I'm afraid this didn't work, either. |
|
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 GitHub<https://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2049784150__;Iw!!Mak6IKo!Je2_4CipPQ3gp8VppUspO3eXypfg4-xZ4BmImfS9nRfYttBjl00rZ7n7rvj4IOxnGdJjgLRhGcHYPSX_pTqwe2RF7Y3q-A$>, or unsubscribe<https://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: ***@***.******@***.***>>
|
|
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. |
|
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). |
|
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 GitHub<https://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2050512619__;Iw!!Mak6IKo!JT_3oWMgZAALm3wW54-RqVHr_JylazEAGdg3bHGkXEhq3mxNJiVeUfvpvvR2iEeGcyquRZL9U5S1nNExlIR8ktL9i0A0iA$>, or unsubscribe<https://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: ***@***.******@***.***>>
|
|
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 WARNING : Resampling found -129405.02% of the output grid 'mesoscale_m1' covered. Will skip producing this product: C03 |
|
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: This maps to the following X/Y coordinates in the geostationary projection: 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... |
|
@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 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 |
|
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)
@mraspaud<https://urldefense.com/v3/__https:/github.com/mraspaud__;!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI6-ktrUDeQ$> @pnuu<https://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 GitHub<https://urldefense.com/v3/__https:/github.com/ssec/polar2grid/issues/691*issuecomment-2065585772__;Iw!!Mak6IKo!LzrA9kJEnDyNXfwLH8c_hcgbESsSduMyXmXlbdY_qqfmIDc8osLn6DLcpvWHXmZcxICyIEdimESM1mflyD1gI6-17Qf4IA$>, or unsubscribe<https://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: ***@***.******@***.***>>
|
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. |
|
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? |
|
@mraspaud You are not missing anything, but @pnuu is right too that
|
|
Oh, right. I missed that the other use is for lonlats. |
|
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 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 |
|
Your last solution looks like a good compromise. |
|
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 |
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:
The text was updated successfully, but these errors were encountered: