Comparing odc.stac.load and stackstac for raster composite workflow

Data
22

Just joining this thread, but I think this point is incredibly important. You need to know how coords are specified in the metadata! Nice job illustrating the fuzziness in the example!

23

I ran a similar test for stackstac.stac and it appears that this difference does not apply to xy_coords?

image
image684×878 96.7 KB

The four images above were built from code like this:

l30_stack = ss.stack(itemsL30, 
                    assets = ['B05'], 
                    epsg = epsg, # epsg (3979 if set to poly) # epsg (32614 if set to UTM)
                    resolution = 30, 
                    bounds = bbox, 
                    resampling = Resampling.bilinear, # bilinear # cubic # nearest
                    chunksize = (1, -1, -1, -1), 
                    xy_coords = 'topleft', # center # topleft
                    gdal_env = gdalEnv)

With just changes to xy_coords and epsg (keeping to local UTM or reprojecting to national proj, as I do in my work). All of these used bilinear resampling, but same patterns held true for cubic and nearest neighbor.

As you can see, xy_coords does not impact “blurriness” like anchor does. The only change is the values are shifted half a pixel between topleft and center. The blurriness in this case occurs when reprojecting from the local UTM grid to a national projection - as would be expected (and required for wide-area analysis when you need everything projected onto the same pixel grid).

So that begs the question, is topleft or center the “correct” pixel placement? I though it would be topleft for HLS, as noted in the User Guide and mentioned by @Andrew_Mullen. But it actually appears that center is most accurate to the real world:

image
image872×794 62 KB

In these comparisons, left is topleft and right is center. These are using the national projection, not UTM, but the patterns are the same.

So I guess the takeaway is that for stackstac that xy_coords = center is the proper way, at least for HLS. I guess this is because it integrates with xarray, which likes pixel center coordinates. And that it is not work the same as anchor does with odc.stac.load?

EDIT: I also played around some with snap_bounds = False in stackstac.stack and it does impact pixel values, but I can’t really tell if its preferred over the default (snap_bounds = True) if already reprojecting?

1 Like