Skip to content

Instantly share code, notes, and snippets.

@mdsumner

mdsumner/netcdf_to_grid.md

Created December 12, 2024 06:49
Show Gist options
  • Save mdsumner/aff0bc4155535f40dfb85d896e5d6420 to your computer and use it in GitHub Desktop.
Save mdsumner/aff0bc4155535f40dfb85d896e5d6420 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
netcdf_to_grid.md

From here https://planetarycomputer.microsoft.com/dataset/sentinel-5p-l2-netcdf#Example-Notebook

we get

import cartopy.crs as ccrs
import fsspec
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import planetary_computer
import pystac_client
import xarray as xr

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

longitude = 79.109
latitude = 22.746

geometry = {
    "type": "Point",
    "coordinates": [longitude, latitude],
}

search = catalog.search(
    collections="sentinel-5p-l2-netcdf",
    intersects=geometry,
    datetime="2023-04-02/2023-04-03",
    query={"s5p:processing_mode": {"eq": "OFFL"}, "s5p:product_name": {"eq": "hcho"}},
)
items = list(search.get_items())

items[0].assets["hcho"].href
'https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D'

In linux I need this option to enable NETCDF (rather than HDF5 driver)

docker run --rm -ti  --security-opt seccomp=unconfined   ghcr.io/mdsumner/gdal-builds:rocker-gdal-dev-python  bash

which gives a massive list of subdatasets

 gdalinfo -if NetCDF "/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D"

I can't find a nice way to organize that output from in R, so I've just picked one

dsn <- "NETCDF:\"/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D\":/PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_temperature"

## gdalinfo on that shows that Geolocation Arrays are present
#Geolocation:
#  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
#  X_DATASET=NETCDF:"/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D":/PRODUCT/longitude
#  X_BAND=1
#  Y_DATASET=NETCDF:"/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D":/PRODUCT/latitude


library(gdalraster)
## use the auto-target grid heuristics
warp(dsn, "/vsimem/surface_temperature.tif", t_srs = "EPSG:4326")

terra::rast("/vsimem/surface_temperature.tif")

class       : SpatRaster
dimensions  : 4578, 9160, 1  (nrow, ncol, nlyr)
resolution  : 0.03930033, 0.03929175  (x, y)
extent      : -180, 179.9911, -89.92088, 89.95674  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source      : surface_temperature.tif
name        : surface_temperature

Oh, so now I realize these are nearly-global datasets. So use this

library(terra)

dsn <- "NETCDF:\"/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D\":/PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_temperature"

project(rast(dsn), rast(res = 0.04), by_util = TRUE)

class       : SpatRaster
dimensions  : 4500, 9000, 1  (nrow, ncol, nlyr)
resolution  : 0.04, 0.04  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
source      : spat_1e14a27bd4e_481_o8up7yuDdJ6XHBX.tif
name        : surface_temperature
min value   :            204.6796
max value   :            310.6992
time        : 2023-04-03 UTC
@mdsumner
Copy link
Author

Use literally any target raster there, so we can do

library(terra)

dsn <- "NETCDF:\"/vsicurl/https://sentinel5euwest.blob.core.windows.net/sentinel-5p/TROPOMI/L2__HCHO__/2023/04/03/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423/S5P_OFFL_L2__HCHO___20230403T063304_20230403T081434_28345_03_020401_20230404T225423.nc?st=2024-12-11T06%3A19%3A34Z&se=2024-12-12T07%3A04%3A34Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-12-07T09%3A00%3A14Z&ske=2024-12-14T09%3A00%3A14Z&sks=b&skv=2024-05-04&sig=tJJZa0pq23QgpI2d0Do33kTNK10ELlUXEQlv24JyCfg%3D\":/PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_temperature"

r <-  project(rast(dsn), rast(ext(-1, 1, -1, 1) * 1e7, crs = "+proj=laea +lon_0=70", res = 25000), by_util = TRUE)
plot(r)

image

happy to explore more, this is just what I found. I'm using very recent gdal but we can probably arrange something for earlier versions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment