Conversion from netcdf to raster (tif) changes resolution and extent?

user_13467482 0
Viewed 279

I'm trying to convert a netcdf file into raster (tif) format. I have created a script and it worked fine a while ago. But now, when I try to use the same simple script with a different file, the resolution changes from 0.5 x 0.5 to 0.5 x 0.5263158. Also the extent moves from:

-100.25, -73.25, 28.75, 48.75

to

-100.5, -73, 28.48684, 49.01316

I have also tried with different raster packages in R, yet they come back with a message saying that cells are not equally spaced. It could be very well an issue with the file (attached here), but I could not see where and how.

Code for reproduction:


# load netcdf file
import xarray as xr 
import rioxarray

xds = xr.open_dataset('output_shocks_us/hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044.nc')
xds = xds.rename({'lat':'y','lon':'x', 'time':'band'})

# Add CRS
xds.rio.write_crs("epsg:4326", inplace=True)

# Convert to geotiff
xds["yield-soy-noirr"].rio.to_raster('hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044_test.tif')
rio = xr.open_rasterio("hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044_test.tif")

print(xds)
print(rio)

The full results are:

print(xds)
<xarray.Dataset>
Dimensions:          (y: 39, x: 55)
Coordinates:
    band             int64 2025
  * y                (y) float64 28.75 30.25 30.75 31.25 ... 47.75 48.25 48.75
  * x                (x) float64 -100.2 -99.75 -99.25 ... -74.25 -73.75 -73.25
    spatial_ref      int32 0
Data variables:
    yield-soy-noirr  (y, x) float64 nan nan nan nan nan ... nan nan nan nan nan
Attributes:
    grid_mapping:  spatial_ref

############
print(rio)
<xarray.DataArray (band: 1, y: 39, x: 55)>
array([[[     nan,      nan, ...,      nan,      nan],
        [     nan,      nan, ...,      nan,      nan],
        ...,
        [     nan, 0.672842, ...,      nan,      nan],
        [     nan,      nan, ...,      nan,      nan]]])
Coordinates:
  * band     (band) int32 1
  * y        (y) float64 28.75 29.28 29.8 30.33 30.86 ... 47.17 47.7 48.22 48.75
  * x        (x) float64 -100.2 -99.75 -99.25 -98.75 ... -74.25 -73.75 -73.25
Attributes:
    transform:      (0.5, 0.0, -100.5, 0.0, 0.5263157894736842, 28.4868421052...
    crs:            +init=epsg:4326
    res:            (0.5, -0.5263157894736842)
    is_tiled:       0
    nodatavals:     (nan,)
    scales:         (1.0,)
    offsets:        (0.0,)
    descriptions:   ('yield-soy-noirr',)
    AREA_OR_POINT:  Area
    grid_mapping:   spatial_ref
terra python-xarray netcdf r python
1 Answers

You're probably not waiting for anyone to answer here anymore, but I'd like to provide an attempt of an answer hereby:

When reading the netCDF file using terra, several problems are evident. There is no coordinate reference system provided, no extent defined and your resolution differs in direction of x and y. Like you said in your question, just wanted to confirm.

nc <- terra::rast("hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc")
nc
#> Error in R_nc4_open: Invalid argument
#> Warning: [rast] GDAL did not find an extent. Cells not equally spaced?
#> class       : SpatRaster 
#> dimensions  : 39, 55, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01818182, 0.02564103  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc:yield-soy-noirr 
#> varname     : yield-soy-noirr 
#> name        : yield-soy-noirr

terra throws a warning, not an error, and nevertheless creates a SpatRaster object.

Since we know the dimensions of your raster (x: 55, y: 39) and the resolution (0.5°), I suspect - please prove me wrong! - that this is inconsistent with the extent provided assuming the lower left corner of your bounding box is correct:

# x
seq(from = -100.25, by = 0.5, length.out = 55 + 1) |> range()
#> [1] -100.25  -72.75

# y
seq(from = 28.75, by = 0.5, length.out = 39 + 1) |> range()
#> [1] 28.75 48.25

... and the result after manually adjusting attributes looks quite promising since your resolution now has the expected values:

terra::ext(nc) <- c(-100.25, -72.75, 28.75, 48.25)
terra::crs(nc) <- "epsg:4326"

nc
#> class       : SpatRaster 
#> dimensions  : 39, 55, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -100.25, -72.75, 28.75, 48.25  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc:yield-soy-noirr 
#> varname     : yield-soy-noirr 
#> name        : yield-soy-noirr

The resulting SpatRaster object can be written to disk using:

terra::writeRaster(nc, "hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.tif")
user_117092960
Related