Copernicus 30m DEM support

[uchchwhash]

Hi @mwheeler.

I just created a new branch for you to play with: https://github.com/OpenDataCubePipelines/ard-pipeline/tree/cop-dem-stash. Go nuts!

Essentially, the task here is to develop code that lets the ARD pipeline use the Copernicus 30m DEM as is. The DEM we currently use is based on SRTM and is in a mozaic. Whereas the CopDEM comes in 1deg-by-1deg tiles. I think you would want to develop it to be able to download from the s3://copernicus-dem-30m/ bucket. Would you please keep the download code separate from the processing code maybe? Not critical.

You can create a stand-alone python script to do this if you want. Integration into our codebase is not necessary.

The relevant code for DEM reading (it says it supports CopDEM, but that's really a mozaic that we created) https://github.com/OpenDataCubePipelines/ard-pipeline/blob/cop-dem-stash/wagl/dsm.py#L45

This function takes an acquisition object, you can create it like:

from wagl.acquisition import acquisitions

acquisition = acquisitions(PATH_TO_YOUR_LEVEL1_DATASET)

With this, you should be able to use much of the code in get_dsm. Roughly what it does is:

• create a geobox (an object that represents a bounded grid (like an image) with origin, pixelsize, shape, and CRS info)
• expand the geobox by some specified buffer (to be able to calculate shadows cast by surrounding elevation)
• read the data into an array with the same geospatial extent as the expanded geobox (this is where we need to re-engineer things)
• re-project into the target CRS and re-sample (I think you can just reuse the code there)
• apply a gaussian filter (not super important, and you can reuse the code)
• write to an hdf5 group (you can write to a tiff, and since you already have the target geobox (dem_geobox in the code I think), you can draw inspirations from https://github.com/OpenDataCubePipelines/ard-pipeline/blob/cop-dem-stash/wagl/data.py#L120 if you want)

So, the key job here is the read the data into an array. We would not have the luxury of calculating one geobox in the DEM CRS to read the data from anymore, but fortunately read_subset only needs extends in the DEM CRS (angular coordinates in our case). I think you can experiment with mozaicing first and then read a subset, or reading overlapping subsets from each tile and then mozaicing them, whichever is convenient. It is not necessary to use read_subset but it supports reading from rasterio datasets if that helps.

In the code, subs is the array of DEM data in DEM coordinates, which we reproject and resample immediately after to the target coordinates, which is called dsm_data. Note that the desired output is really dsm_data. Mix and match the data read part and the reproject/resample part as you see fit.

And please feel free to reach out if you need any clarifications.

[uchchwhash]

Just clarifying something that I just got bitten by:

from wagl.acquisition import acquisitions

# the acquisitions function actually returns an `AcquisitionContainer` object, not an `Acquisition` object
container = acquisitions(PATH_TO_YOUR_LEVEL1_DATASET)

# this is what you want to pass to `get_dsm`
acquisition = container.get_all_acquisitions()[0]

or something like that.

[mwheeler]

Just clarifying something that I just got bitten by:

from wagl.acquisition import acquisitions

# the acquisitions function actually returns an `AcquisitionContainer` object, not an `Acquisition` object
container = acquisitions(PATH_TO_YOUR_LEVEL1_DATASET)

# this is what you want to pass to `get_dsm`
acquisition = container.get_all_acquisitions()[0]

or something like that.

Yep, all good 👍 I haven't had any issues figuring out the wagl functions so far / there's already a prototype of the DEM acquisition here: https://github.com/OpenDataCubePipelines/ard-pipeline/commit/3929a79ad90b2ec4d377bd40d687769d656f4f36

I still have a weird bug to solve (west most column of tiles don't seem to get mosaiced, unsure why yet - need to manually go through the process to compare if maybe it actually is correct and visually it just doesn't meet what I expected with a 1deg border, or if it's actually wrong (it looks wrong... way too much space on eastern side, not enough on western side - despite north/south border looking fine with same math...)).

I've put that to the side for a day though & I'm already working on the other ECMWF ancillary data (should be done sometime this week depending on my other workload).

Depending on any bugs I have to fix, I'm hoping I can give you a nice polished/documented/tested bunch of code next week 🤞 (for all 4 ancillary, DEM+watervapour+ozone+aerosol)

[uchchwhash]

Thanks @mwheeler , closing this.