How do I open one of the goes16 netcdf files using rasterio?

[john-a-m]

I'm able to open the netcdf files with xarray, however xarray is slow when opening the dataset and when doing reprojections. I'm trying out rasterio to see if it's faster. Here is the code I've come up with so far.

import rasterio
from affine import Affine

transform = Affine(0.02906874414461869, 0.0, -124.74358263413838, 0.0, -0.029068744144618696, 49.39151116004554)
path = "netcdf:///root/summary24hr/noaa-goes16/ABI-L2-RRQPEF/2024/124/22/OR_ABI-L2-RRQPEF-M6_G16_s20241242200204_e20241242209512_c20241242210000.nc:RRQPE"

with rasterio.open(path, 'r+') as f:
    #f.write_transform(transform)
    f.transform = transform

When opening the dataset I get a NotGeoreferencedWarning. When I try and set the transform I get an error

Traceback (most recent call last):
  File "/root/rasteriotest.py", line 15, in <module>
    f.transform = transform
    ^^^^^^^^^^^
  File "rasterio/_base.pyx", line 721, in rasterio._base.DatasetBase.transform.__set__
  File "rasterio/_io.pyx", line 1655, in rasterio._io.DatasetWriterBase.write_transform
TypeError: not all arguments converted during string formatting

Here is a link to the file I am working with https://noaa-goes16.s3.amazonaws.com/ABI-L2-RRQPEF/2024/124/22/OR_ABI-L2-RRQPEF-M6_G16_s20241242200204_e20241242209512_c20241242210000.nc

[snowman2]

Writing netcdf (or other multidimensional) files is not necessarily something I would recommend doing with rasterio. Reading works, but I would avoid writing if possible.

Since you mentioned xarray, I am guessing that you have seen this already: https://gis.stackexchange.com/questions/345650/how-do-i-add-projection-to-this-netcdf-file-satellite

[john-a-m]

I discovered a third way (https://gis.stackexchange.com/a/414528) that gives yet another shape (1030, 2444).

coord_ref_sys = CRS.from_cf(dataset.goes_imager_projection.attrs)
dataset.rio.write_crs(coord_ref_sys.to_string(), inplace=True)
bbox = box(-124.725839, 24.498131, -66.949895, 49.384358)
feature_proj = rasterio.warp.transform_geom(
    CRS.from_epsg(4326),
    coord_ref_sys,
    bbox
)
print(feature_proj)
ds = dataset.rio.clip([feature_proj])
numpy_array = ds[["RRQPE"]].to_array().to_numpy()[0]
print(numpy_array.shape)

[snowman2]

This should do the same thing in a simpler way:

coord_ref_sys = CRS.from_cf(dataset.goes_imager_projection.attrs)
dataset.rio.write_crs(coord_ref_sys, inplace=True)
bbox = box(-124.725839, 24.498131, -66.949895, 49.384358)
dataset = dataset.rio.clip([bbox], crs="EPSG:4326")
print(dataset.rio.shape)

[snowman2]

I would imagine that this method is the most efficient as it reduces the amount of data you reproject:

dataset = dataset.rio.clip_box(
    minx=-124.725839,
    miny=24.498131,
    maxx=-66.949895,
    maxy=49.384358,
    crs="EPSG:4326"
)
dataset = dataset.rio.reproject("EPSG:4326")
dataset = dataset.rio.clip_box(
    minx=-124.725839,
    miny=24.498131,
    maxx=-66.949895,
    maxy=49.384358,
    crs="EPSG:4326"
)
print(dataset.rio.shape)

The reason for doing the clip twice is that the first clip transforms the boundary to clip into the raster projection. This often results in an expansion of the boundary.

See: https://corteva.github.io/rioxarray/stable/examples/transform_bounds.html

Performing the clip on the raster post reprojection will use the boundary in its original format without expansion.

[john-a-m]

Something like what you've written is what I'll end up going with. I think I can just drop the reproject all together. Then when I draw the raster to leaflet I'll just compute the pixels manually xWidth = ((xmax - xmin) / array_width) and then xOffset = xIndex * xWidth. I'm still somewhat inexperienced with the GIS stuff, so I always try to verify the outputs of my programs with either gdal or qgis. To that end I've asked a question https://gis.stackexchange.com/questions/480867 if you're interested. In any case thanks for all your help and patience.

[john-a-m]

I wrote my first reply while you were writing about clipping twice. I'll give that a try.