Mosaic / Reproject methods

[bbbarnes]

I'm trying to decipher what happens under the hood of the GPT Mosaic and Reproject routines (SeaDAS 7.4).

Here's what I've done that's got me to the point of totally confused: Start with A2016316180000.L1A_LAC -> modis_GEO.py -> modis_L1A_extract.py (24 to 26 Lat, -83 to -80 Lon) -> modis_L1B.py -> l2gen with all defaults except resolution=250. Note this extraction is very near the scan edge to highlight my problems. Also, I'm showing reprojections at 250m resolution, even though I recognize that even Rrs(645) does not have that fine resolution near the scan edge. The same effects are seen for HKM resolution projections, they're just not as obvious...

Anyway, I can mosaic this file to the default projection (with the specified bounds above and 0.0022727 pixel resolution [approx 250m] - see attached .xml), and output Rrs(645) to .png [mosaic_nearest.png]. If I retry this after changing the resampling method in the .xml to either 'Bicubic' or 'Bilinear', I get the exact same image [mosaic_bilinear.png; mosaic_bicubic.png]. I don't see how this is possible, unless the bilinear / bicubic methods are parameterized such that they do not "find" any neighbors for interpolation around each pixel in the output grid. Or perhaps under certain conditions that resampling option is ignored?

Interestingly, I tried the same thing with Reproject (using the GUI, in case my mosaic .xml was messing things up). The Nearest reproject [reproj_nearest.png] looked almost exactly the same as the Nearest mosaic, with just a few speckles filled. The Bilinear and Bicubic options in reproject produced exactly the same values as the nearest for water pixels, but the land and cloud NaNs were dilated a bit [reproj_bilinear.png, reproj_bicubic.png].

Also, per other forum posts, I reran the Nearest reproject after attaching the GEO coding [reproj_nearest_geocode_attached.png]. Here, I can see some changes in the data as well as the masking.

On top of all of this, I've tried several other softwares to reproject these data (using just the latitude, longitude, and Rrs_645 datasets in the L2) to a cylindrical equidistant projection using nearest neighbor and always get the same answer, which is quite different from everything above [pyresample.png]. Note that there are many apparent errors, primarily around the land/cloud edges. To validate that this method was actually nearest neighbor resampling, I tried was physically calculating the nearest neighbor from the lat/lon for each grid point in the output array, and I get the same answer as pyresample.

So, my questions are: What are the (likely fixed) input parameters for the Nearest, Bilinear, and Bicubic methods (i.e., radius of influence, number of neighbors to consider for Bilinear and Bicubic, etc)? Why do these look the same regardless of the selected sampling method (or essentially the same, for reproject)? Most importantly, what latitude and longitude grids are used in mosaic/reproject, since it does not appear to be the latitude and longitude datasets in the L2?

Any documentation / help on this would be greatly appreciated. Thanks in advance.
-brian
attachment 1

[gnwiii]

It would be interesting to know if SNAP 6 (pre-release test version) behaves differently.  Whenever SeaDAS 7 gpt gets annoying I try the SNAP version.  As for documentation, GPT is from ESA and Brockmann so look for Mosaic isses there.   Also, the sources are available.  I find SNAP documentation is generally better than the BEAM versions.  Some warp problems in gdal are due to lack of precision with 32-bit floats -- SNAP went to 64-bit to support higher resolutions.

I've been down this rabbit hole before -- some of the geo-processing code in BEAM and SeaDAS 7.4 came from GeoTools, which is heavily influenced by GDAL, so you could see how gdalwarp handles your files.  In the past, some "interesting" behaviours of BEAM could be traced back to known bugs in GeoTools.  With SNAP,  I'm not sure if GeoTools are still used, but that may not matter, as GeoTools'  warp operation wraps the JAI implementation.

[OB SeaDAS - knowles]

Basically, regarding the source file, bilinear uses a grid size of 2x2 and bicubic uses a grid size of 4x4.
(see link or also available within the SeaDAS internal help): https://seadas.gsfc.nasa.gov/help/general/ResamplingMethods.html

The reason you are not seeing a difference between these methods is that the resolution of the source is much lower than that of the output mapped image due to the high sensor zenith angle in your specific region for this specific file.  The choice of nearest, bilinear, and bicubic only makes a difference when the source resolution is comparable or higher to that of your output mapped resolution file.

To clarify this scan angle dependence on actual ground resolution (some users may not be aware of this issue):
If you look at the attached image which shows the bounds of your level2 file (in white) and your mapped region (in red) you will see that the region is to the westernmost portion of the level2 file.  This is at the end of scan (Aqua scans east to west, but scan direction doesn't matter for this issue).  Both MODIS instruments scan from -55 to +55 degrees (with respect to the instrument) with translates to approx -65 to +65 degree sensor zenith angle (in other words the ground angle).  When referring to resolution in the satellite native pixel resolution raster files (level-0, level1 and level2), the "name" of the resolution is in general that of the (center) nadir (0 degrees) sensor zenith pixel.  However, at the ends of scan (the sensor zenith angle is very high), the actual footprint of the detector's image on the Earth is, in the case of MODIS, about 10 times the size (in area) as that of a pixel at nadir.  In your specific case, the pixels of your source file in your projected region are approximately 500m x 1250m (the footprint stretches longer west-east than north-south).  So you are actually projecting to a much higher resolution of 250m and consequently little difference is observed between the choice of resampling methods. 

Danny

[bbbarnes]

Thanks Danny. Given a bilinear grid size of 2x2, I'm not that surprised that there's no difference between the Nearest and Bilinear. However, I would expect some difference from Bicubic (4x4 grid) - even for a footprint at the scan edge (a 4 x 4 grid would cover 1km x 1km, which is bigger than the footprint size). The outputs are exactly the same (difference = 0), so it's not just that there is "little difference".

The bigger question for me, though, is what are the input lat / lon fields used for the warping? True nearest neighbor (physically calculating the closest lon/lat for each output grid pixel) using the L2 latitude and longitude datasets yields the pyresample.png image above, which looks very dissimilar to the GPT nearest.

@gnwiii- I've tried gdalwarp in the past without much success. By default, it warps just using the boundary coordinates from the L2, so the scan line banding is retained [gdalwarp.png]. Using the -geoloc option (which should force it to use the included lat/lon datasets) throws an error. You can get around the error by converting everything to .vrt files and manually changing the geoloc parameters in them, then run with -geoloc, but the output leaves something to be desired [gdalwarp_geoloc.png]. I'll give SNAP a shot next I guess.

-brian

[gnwiii]

See How to georeference MODIS Terra SST? for messy example using GDAL with level-2 lon/lat fields extracted to a text file and the thin plate splines option of gdalwarp.  GDAL has been improving, so if your current version gives you grief it is worth the effort to a) build a current version and b) find a way for the new version to coexist with the old version so you don't break things linked against the old version.

[OB SeaDAS - knowles]

Regarding " (a 4 x 4 grid would cover 1km x 1km, which is bigger than the footprint size)"
The grid used is that of the source file not the target "reprojected" file.   When you reproject to 250m, the adjacent pixels of the source file are orders of magnitude away from the nearest neighbor pixel of the source file pixel and consequently should have negligible weighting impact on the interpolation regardless of chosen method.  Perhaps the difference will show up way down in the decimal places?

Also since you are using the default projection (geo lat lon WGS84) the projected pixel sizes are not equal area so your output file pixels are not 250m (varies with latitude) but certainly smaller than the source pixels.

I'm not sure what you mean by warping.  The end of scan pixels certainly have a warped ground projection as well as significantly overlap eachother.  There is also a "bow-tie" effect at these large scan angles between scans causing the striping effect (where the latitude doesn't continually increase/decrease across the vertical path through the raster pixels of the level2 file).  The SeaDAS reader makes corrections for this effect.  Attaching geocoding skips this correction.  Perhaps that may be the source of some of the discrepancies you are seeing.

Danny

[bbbarnes]

Thanks Danny.

I tried at ~8km resolution - same xml as above, but using:
         <pixelSizeX>0.072727272</pixelSizeX>
         <pixelSizeY>0.072727272</pixelSizeY>

Again, the outputs of Nearest, Bilinear, and Bicubic are identical (at least to floating point precision, which is what Mosaic defaults to).  Same thing with HKM, 1KM, 2KM and 4KM...

My apologies about the term "warping" (I had gdalwarp on my mind). I meant "reprojection".

Back to the QKM projection - my problem is that, for a given output grid pixel (center lat/lon), the "nearest neighbor" in the L2 as reported by mosaic is not always the closest input latitude and longitude combination from the lat / lon datasets in the L2. I know the pyresample image above shows this, as I've checked it by manually calculating the distances on a per-pixel basis. Perhaps the difference is a result of the bowtie correction, but I cannot replicate the pyresample behavior by attaching the geo-coding (even varying the search radius).

So, what does the reader feed to GPT? Are some lat/lons different from the L2 datasets? Are some lat/lons ignored? Are only ground control points passed?

Thank you for being so attentive.
-brian

[OB SeaDAS - knowles]

Brian,

In order to better help you with this and separate out the end of scan issues from nadir issues you might want to try the same file but using a region at nadir (in the  center file), say (21N 19S -71W -68E).  Let us know if you see the same type of discrepancies between your calculated nearest neighbor pixel and the one Reproject chooses.  We always appreciate user feedback, integrity testing, kicking the tires, etc. :)

Of note is that the Mosaic Tool actually internally calls the Reproject GPT operator so testing Reproject is sufficient.

The logic of the bow-tie correction is likely different than what you might be applying in your nearest pixel calculations.  In a nut shell, if a file invokes BowtiePixelGeoCoding, the projection transform first finds the nearest scan to the target lat/lon ___location (by comparing center of scan lat/lon to target lat/lon).  Then it chooses the nearest pixel (within this chosen scan) to the target lat/lon ___location.  MODIS scan widths are (LAC=10 detectors/pixels; HKM=20 detectors/pixels; QKM=40 detectors/pixels).

Some specific answers to your questions:
What does the reader feed to GPT? Essentially the knowledge of whether the source file contains a known CRS (coordinate reference system), or if its unmapped.  If unmapped then whether to use the bowtie correction.
Are some lat/lons different from the L2 datasets? No
Are some lat/lons ignored? Yes, per the bowtie logic, regarding oversampled areas, the pixel of the "better" scan is used and the one of the "worse" scan is ignored. 

Regarding your observation "the outputs of Nearest, Bilinear, and Bicubic are identical", yes I do believe there is (currently) logic which forces "Nearest Neighbor" resampling to be used under certain conditions, one being that the file is unmapped (has no CRS).  We will need to take a look into this in order to positively confirm this in indeed the behavior.

Danny

[bbbarnes]

Thanks again Danny. That's just what I needed to know!

If I duplicate the bowtie logic, the answer I get from pyresample is the same as the GPT output for the area near the scan edge. Yay!

Note I tried the GPT/pyresample comparison for the region you proposed at the center of the image. The pyresample image looked correct (didn't show any of the discontinuities I was seeing for the scan edge region), but I couldn't compare directly to the GPT output because my attempt to run Reproject threw our old friend, the linestring error...

Either way, I think you've addressed my concerns about what exactly GPT is doing for Nearest. Still not clear on why Nearest, Bilinear, and Bicubic are the same, but I'm assuming that GPT is reverting to default Nearest for some reason.

Thanks again for all of your help.
-brian