#!/usr/bin/env python3 # Convert decadal NetCDF variables on a curvilinear grid to multiband # EPSG:3031 Cloud Optimised GeoTIFFs. # # Install: # pip install numpy # # GDAL Python bindings must also be installed and compatible with the # installed GDAL binaries (on Windows, usually via a matching wheel). from pathlib import Path import numpy as np import xarray as xr from osgeo import gdal gdal.UseExceptions() # ================================================== # CONFIG # ================================================== BASE_DIR = Path.cwd() CONFIG = [ { "nc": "ecoindex_historical_future.nc", "vars": [ ( "ecoindex_hemispheric_decadal", "AEV_hemispheric_decadal_multiband.tif", ), ( "ecoindex_regional_decadal", "AEV_regional_decadal_multiband.tif", ), ], }, ] # EPSG:3031 output grid. XMIN, YMIN, XMAX, YMAX = -4500000, -4500000, 4500000, 4500000 RES = 25000 # 25 km NODATA = -9999 COG_OPTIONS = [ "COMPRESS=DEFLATE", "BLOCKSIZE=512", "BIGTIFF=IF_SAFER", ] # ================================================== # FUNCTIONS # ================================================== def make_geolocation_files(nc_path: Path): """ Extract TLAT and create a corrected TLONG raster for GDAL geolocation. """ stem = nc_path.stem tlat_path = BASE_DIR / f"TLAT_raw_{stem}.tif" tlong_path = BASE_DIR / f"TLONG_fixed_{stem}.tif" if not tlat_path.exists(): gdal.Translate( str(tlat_path), f'NETCDF:"{nc_path}":TLAT', format="GTiff", creationOptions=[ "TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512", "COMPRESS=DEFLATE", ], ) if not tlong_path.exists(): src = gdal.Open(f'NETCDF:"{nc_path}":TLONG') if src is None: raise RuntimeError(f"Could not read TLONG from {nc_path}") array = src.ReadAsArray().astype(np.float32) array[~np.isfinite(array)] = np.nan array[array < -180.0] += 360.0 out = gdal.GetDriverByName("GTiff").Create( str(tlong_path), src.RasterXSize, src.RasterYSize, 1, gdal.GDT_Float32, options=[ "TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512", "COMPRESS=DEFLATE", ], ) out.SetGeoTransform((0, 1, 0, 0, 0, 1)) band = out.GetRasterBand(1) band.WriteArray(array) band.SetNoDataValue(np.nan) out = None src = None return tlat_path, tlong_path def build_slice_vrt( nc_path: Path, variable: str, band_number: int, tlat_path: Path, tlong_path: Path, vrt_path: Path, ): """ Create a single-band VRT for one time slice and attach geolocation arrays. """ gdal.Translate( str(vrt_path), f'NETCDF:"{nc_path}":{variable}', format="VRT", bandList=[band_number], ) ds = gdal.Open(str(vrt_path), gdal.GA_Update) ds.SetMetadata( { "X_DATASET": tlong_path.name, "X_BAND": "1", "Y_DATASET": tlat_path.name, "Y_BAND": "1", "PIXEL_OFFSET": "0", "PIXEL_STEP": "1", "LINE_OFFSET": "0", "LINE_STEP": "1", }, "GEOLOCATION", ) ds = None def warp_slice(vrt_path: Path, output_path: Path): """ Reproject one time slice to the fixed EPSG:3031 grid. """ gdal.Warp( str(output_path), str(vrt_path), options=gdal.WarpOptions( format="GTiff", srcSRS="EPSG:4326", dstSRS="EPSG:3031", geoloc=True, xRes=RES, yRes=RES, outputBounds=[XMIN, YMIN, XMAX, YMAX], resampleAlg="near", dstNodata=NODATA, creationOptions=[ "TILED=YES", "COMPRESS=DEFLATE", "BIGTIFF=IF_SAFER", ], ), ) def stack_to_cog( band_files: list[Path], band_labels: list[str], output_path: Path, ): """ Stack projected time slices into a multiband COG and set band descriptions. """ vrt_path = output_path.with_suffix(".vrt") gdal.BuildVRT( str(vrt_path), [str(path) for path in band_files], separate=True, ) ds = gdal.Open(str(vrt_path), gdal.GA_Update) ds.SetMetadata({}) for band_number, label in enumerate(band_labels, start=1): band = ds.GetRasterBand(band_number) band.SetMetadata({}) band.SetDescription(label) ds = None gdal.Translate( str(output_path), str(vrt_path), format="COG", creationOptions=COG_OPTIONS, ) vrt_path.unlink() def get_time_values(variable, dataset): """ Find the time/decade coordinate associated with a variable. """ for name in variable.coords: lower = name.lower() if any(term in lower for term in ["time", "year", "decade"]): return variable.coords[name].values for dimension in variable.dims: lower = dimension.lower() if all(term not in lower for term in ["lat", "lon", "x", "y"]): if dimension in dataset.coords: return dataset.coords[dimension].values if dimension in dataset: return dataset[dimension].values return np.arange(variable.shape[0]) def make_decade_label(value, index): """ Convert a time coordinate to a decade label such as 2050-2060. """ try: decade_start = int(value) except (TypeError, ValueError): try: year = int(str(value)[:4]) decade_start = (year // 10) * 10 except (TypeError, ValueError): decade_start = index * 10 return decade_start, f"{decade_start}-{decade_start + 10}" # ================================================== # MAIN # ================================================== def main(): for config in CONFIG: nc_path = BASE_DIR / config["nc"] if not nc_path.exists(): raise FileNotFoundError(nc_path) tlat_path, tlong_path = make_geolocation_files(nc_path) with xr.open_dataset(nc_path) as dataset: for variable_name, output_name in config["vars"]: variable = dataset[variable_name] times = get_time_values(variable, dataset) band_files = [] band_labels = [] for index, time_value in enumerate(times): decade_start, label = make_decade_label( time_value, index, ) vrt_path = ( BASE_DIR / f"{variable_name}_{decade_start}.vrt" ) temp_path = ( BASE_DIR / f"{variable_name}_{decade_start}_warp.tif" ) build_slice_vrt( nc_path, variable_name, index + 1, tlat_path, tlong_path, vrt_path, ) warp_slice( vrt_path, temp_path, ) band_files.append(temp_path) band_labels.append(label) vrt_path.unlink() output_path = BASE_DIR / output_name stack_to_cog( band_files, band_labels, output_path, ) for path in band_files: path.unlink() print(f"Created: {output_path.name}") print("Done.") if __name__ == "__main__": main()