3

My code creates a geoTiff from a NetCDF and then I'm trying to compute zonal_stats() from a raster file I created from a NetCDF. I end up getting ValueError: negative dimensions are not allowed at the last step of my short code:

import geopandas as gpd
import rioxarray as rxr
import rasterstats as rs
import xarray as xr
import rioxarray

# Open the world-wide NetCDF
file_nc=xr.open_dataset("wglc_timeseries_05m.nc")
NC_Data = file_nc.rio.write_crs("epsg:4326", inplace=True)

# Slice the first date
export_raster=NC_Data.sel(time=NC_Data["time"].values[1])
export_raster.rio.to_raster("temp.tif") 
reimported_raster = rxr.open_rasterio("temp.tif").squeeze()

# Load african centroids to prepare buffers
centroids_points = gpd.read_file("centroids.json", driver="GeoJSON")
centroids_poly = centroids_points.copy()                
centroids_poly["geometry"] = centroids_points.geometry.buffer(0.1)
centroids_poly.to_file("plot_buffer.shp")

# Collect buffer stats
buffer_light = rs.zonal_stats("plot_buffer.shp",
                                    reimported_raster.values,
                                    affine=reimported_raster.rio.transform(),
                                    stats="count mean")

I have the feeling this might be related to the negative numbers here, but I'm unsure how to fix this.

reimported_raster.rio.transform()
Out[445]: 
Affine(0.08333333333333318, 0.0, -179.99999999999966,
       0.0, 0.0833333333333333, -89.99999999999997)

For replication, here are the downloads of wglc_timeseries_05m.nc and the centroids.json

1 Answer 1

4

Actually, it's related to a positive number, the y pixel dimension should be negative (pixel coords are top->bottom, but map coords are bottom->top), but it isn't. This is why a negative dimension is being reported - i.e. rasterstats thinks the shape is (-2160, 4320) instead of (2160, 4320).

The affine transform isn't being read correctly by the default xarray engine. Using the rasterio engine in following works for me:

file_nc = xr.open_dataset("wglc_timeseries_05m.nc", engine="rasterio")
etc...
reimported_raster.rio.transform()

Affine(0.08333333333333316, 0.0, -179.99999999999966,
   0.0, -0.08333333333333329, 89.99999999999997) # <--- Note fifth element (transform[4]) is negative.

As a bonus, you don't need to write out to a tif on disk if you don't want to:

# Collect buffer stats
buffer_light = rs.zonal_stats("plot_buffer.shp",
    export_raster.to_array().squeeze().values,
    affine=export_raster.rio.transform(),
    stats="count mean")

buffer_light
[{'mean': 0.0, 'count': 5}, {'mean': 0.0, 'count': 4}, {'mean': 0.005268027121201157, 'count': 5}, {'mean': 0.0, 'count': 5}, {'mean': 0.0, 'count': 5}, {'mean': 0.0, 'count': 5} etc...

And with just rioxarray:

# Open the world-wide NetCDF
ds = rxr.open_rasterio("wglc_timeseries_05m.nc")
ds = ds.rio.write_crs("epsg:4326")

# Slice the first date
ds0 = ds.sel(time=ds["time"].values[0])

# Load african centroids to prepare buffers
centroids_points = gpd.read_file("centroids.json", driver="GeoJSON")
centroids_poly = centroids_points.copy()
centroids_poly["geometry"] = centroids_points.geometry.buffer(0.1)
centroids_poly.to_file("plot_buffer.shp")

# Collect buffer stats
buffer_light = rs.zonal_stats(
    "plot_buffer.shp",
    ds0.squeeze().values,
    affine=ds0.rio.transform(),
    stats="count mean")

print(buffer_light)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.