bedmachine5 to geotiff

bedmachine2geotiff.py

import rioxarray as rio
from rasterio.enums import Resampling
import xarray as xr
import numpy as np
ds = rio.open_rasterio('BedMachineGreenland-v5.nc')


props_to_save = ["thickness", "surface", "bed", "errbed"]

ds.thickness.rio.write_nodata(0, inplace=True)
ds.surface.rio.write_nodata(-9999, inplace=True)
ds.bed.rio.write_nodata(-9999, inplace=True)

for prop in props_to_save:
    print(prop)
    ds[prop].rio.to_raster(raster_path=f"bedmachine-v5-{prop}.tif", 
                            driver="COG", compress='LZW', 
                            tiled=True, dtype='int16', windowed=True, predictor=2)

ds.mask.rio.to_raster(raster_path=f"bedmachine-v5-mask.tif", 
                        driver="COG", compress='LZW', 
                        tiled=True, windowed=True, predictor=2)



# afterwards run this to make 1km versions too
# gdal_translate -of GTiff -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "TILED=YES" bedmachine-v5-thickness.tif thickness1km.tif -tr 1000 1000 -r average