example code?

[mdsumner]

do you want some examples? I'm happy to work on several of the cases - rasterio, odc-geo, rioxarray, and maybe the raster_tools and geoutils.

I'm interested to show how I would go about this, and discuss the various idioms available for GDAL (and how those are seen via various downstream packages).

At any rate I'm keen to discuss so I hope this is welcome.

osgeo.gdal

Here I'm using GDAL syntax for "vrt://" which saves us from some more awkward changes on Open() or via VRT otherwise.

## the file is here
#srcurl = "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"

## and it's best to fix the extent, else the warper goes a less efficient pathway (from the imprecise lon,lat arrays)
## https://discourse.pangeo.io/t/fixing-ghrsst-seeking-opinions/3833/6?u=michael_sumner
#dsn = f'vrt:///vsicurl/{srcurl}?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst'

## we need earthdata creds, which I have in a file so
#import os
#os.environ["GDAL_HTTP_HEADER_FILE"] = "/home/mdsumner/earthdata"

## if we don't want to stream, we can ignore everything above (maybe this is in our earthaccess_data/)
dsn = f'vrt://20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst'
## what's our target tile?  (we want native res I expect so we won't supply any other args)
te = [100, -60, 100 + 256/100, -60 + 256/100]  ## probably I need to be more careful to align to source pixel edges

from osgeo import gdal
gdal.UseExceptions()
ds = gdal.Open(dsn)
## we can get a "live" dataset or virtual one "/vsimem" or an actual file
wds = gdal.Warp("/tmp/tile.tif", ds, outputBounds = te)
wds.RasterXSize
#256
wds.RasterYSize
#256
wds.GetDescription()
#'/tmp/tile.tif'
wds.GetGeoTransform()
# (100.0, 0.010000000000000009, 0.0, -57.44, 0.0, -0.010000000000000009)

[mdsumner]

I'm very very new to profiling in python ... but I adapted the above approach to write the single zoom 0 tile from a MURSST to a png. (It won't be fully correct for scaling and so on, but I think it's probably representative of what should be possible).

Interestingly it seems to take more memory to write to PNG compared to TIF which I think is worth exploring.

Write to png (it took much longer to do this from the source url, but it has the same memory footprint).

from memory_profiler import profile

dsn = f'vrt://../earthaccess_data/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst&scale=-26800.000,7988.000&a_srs=EPSG:4326'

from osgeo import gdal
gdal.UseExceptions()

@profile
def run_tile0(dsn): 
  te = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
  ds = gdal.Open(dsn)
  wds = gdal.Warp("/tmp/tile.png", ds, outputBounds = te, dstSRS = "EPSG:3857", width = 256, height = 256, outputType = gdal.GDT_Byte)
  return 1



if __name__ == '__main__':
    run_tile0(dsn)

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    13    129.1 MiB    129.1 MiB           1   @profile
    14                                         def run_tile0(dsn):
    15    129.1 MiB      0.0 MiB           1     te = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
    16    140.9 MiB     11.8 MiB           1     ds = gdal.Open(dsn)
    17   2153.4 MiB   2012.5 MiB           1     wds = gdal.Warp("/tmp/tile.png", ds, outputBounds = te, dstSRS = "EPSG:3857", width = 256, height = 256, outputType = gdal.GDT_Byte)
    18   2153.4 MiB      0.0 MiB           1     return 1

Write to tif uses less memory (not so helpful, but maybe worth writing to /vismem/.tif and then writing that out to PNG).

from memory_profiler import profile

dsn = f'vrt://../earthaccess_data/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst&scale=-26800.000,7988.000&a_srs=EPSG:4326'

from osgeo import gdal
gdal.UseExceptions()

@profile
def run_tile0(dsn): 
  te = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
  ds = gdal.Open(dsn)
  wds = gdal.Warp("/tmp/tile.tif", ds, outputBounds = te, dstSRS = "EPSG:3857", width = 256, height = 256, outputType = gdal.GDT_Byte)
  return 1



if __name__ == '__main__':
    run_tile0(dsn)

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    13    126.8 MiB    126.8 MiB           1   @profile
    14                                         def run_tile0(dsn):
    15    126.8 MiB      0.0 MiB           1     te = [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
    16    138.3 MiB     11.5 MiB           1     ds = gdal.Open(dsn)
    17    572.9 MiB    434.6 MiB           1     wds = gdal.Warp("/tmp/tile.tif", ds, outputBounds = te, dstSRS = "EPSG:3857", width = 256, height = 256, outputType = gdal.GDT_Byte)
    18    572.9 MiB      0.0 MiB           1     return 1

I don't think you should see a 12Gb spike ever.

[maxrjones]

Your contributions are super welcome, thank you @mdsumner! It's helpful to know that fixing the input coordinates can improve efficiency.

I think your results are roughly similar to what I saw using gdal with vsis3. Those preliminary results are in the table at the bottom of https://github.com/developmentseed/warp-resample-profiling/blob/add-prelim-results/02-process-profiling-results.ipynb. The memory usage was manageable, but the wall-time was rather slow. Was that in line with your experiences working with this data?

[mdsumner]

ah I see, I have a lot to do to get across your workings here. The thing that I thought I understood from here https://developmentseed.org/tile-benchmarking/02-run-tests/06-cmr-netcdf.html is that the 12Gb spike in memory (final plot on the page) was the concern.

My general timing benchmark is 3min to convert an entire file to COG (on compute in Australia). To do that one worst zoom=0 tile is about 30-40s, and so I guess that makes sense as each layer is subdivided with many smaller much more manageable tasks we add up to a few minutes.

The timing for that worst tile seems reasonable to me, given that it's from a poorly-structured-for-the-task archival format, 36000x17999->256x256.

[maxrjones]

@mdsumner I haven't been able to get the vsicurl approach from the commented out section of your first comment to work. Would you be willing to take a look at https://github.com/developmentseed/warp-resample-profiling/blob/vsicurl/examples/load-gdal-netcdf-vsicurl.ipynb to let me know where I'm going wrong? I reduced the example to just the failing load component.

Using vsis3 works...I expect I'm doing something wrong with the authentication.

[mdsumner]

Maybe your GDAL version doesn't support those args for vrt:// syntax?

I'll check, and maybe rewrite it so it's generic to older versions.

[mdsumner]

In the second one you aren't setting the extent and crs of the opened dataset. So the warp is still (possibly problematically) reading from the 1D geolocation arrays context. I didn't actually know we could treat this source as s3 so I'll be exploring that also! 👌

[mdsumner]

I saw the note about not being able to use IMERG with rasterio, so here's a way to use multidim mode to cast to classic with permute:

from osgeo import gdal
gdal.UseExceptions()

## local
ds = gdal.OpenEx("3B-HHR-E.MS.MRG.3IMERG.19980101-S000000-E002959.0000.V07B.HDF5", gdal.OF_MULTIDIM_RASTER)
## or remote, for this I set GDAL_HTTP_HEADER_FILE to a file with my "Authorization: Bearer ..."
#ds = gdal.OpenEx("/vsicurl/https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGHHE.07/1998/001/3B-HHR-E.MS.MRG.3IMERG.19980101-S000000-E002959.0000.V07B.HDF5", gdal.OF_MULTIDIM_RASTER)
rg = ds.GetRootGroup()
md = rg.OpenMDArrayFromFullname("/Grid/probabilityLiquidPrecipitation")
## cast to classic dataset, with dimension permute 
##   advantage is we get the grid referencing that the netcdf/hd5 drivers do not (we don't get the crs tho)
cd = md.AsClassicDataset(1, 2)
cd.RasterXSize
#3600
cd.RasterYSize
#1800
cd.GetGeoTransform()
#(-179.99999694739424, 0.09999999830410791, 0.0, -89.99999694654582, 0.0, 0.09999999660727314)
cd.GetSpatialRef()  ## None, sad face


## now we can warp (I'm just plucking an extent from the ether, needs an index or use of outputBoundsSRS for longlat input etc, and I'm setting input crs here)
w = gdal.Warp("out1.tif", cd, outputBounds =  [-1e5, -1e5, 1e5, 1e5], srcSRS = "EPSG:4326", dstSRS = "EPSG:3857", width = 256, height = 0)
w.RasterXSize
#256
w.RasterYSize
#256
w.GetGeoTransform()
#(-100000.0, 781.25, 0.0, 100000.0, 0.0, -781.25)

I'm not sure if the classic api can transpose, maybe it can so then we could have rasterio do that.