Search code examples
pythonrpython-xarrayterrarasterio

Making a empty/template raster using rioxarray / rasterio in Python


I am currently trying to convert R code into Python. I am trying to make a 'template' raster file which I will then "reproject_match" many other rasters to in my workflow so that they all align.

I can do it in R using the terra package, but am unsure how to do so in python using the rioxarray package. I'm hoping someone might be able to help? Here is my R code:

# load lib
library("terra")

# Make an empty raster. By default is has extent of the earth and projection wgs84. Setting the resolution means the rows and cols are automatically calculated.
a <- rast(res = 1)

# Fill it with some random integer values
vals(a) <- as.integer(runif(nrow(a) * ncol(a), 1, 10))

# Write it out to file as an
writeRaster(a, "my_integer_rast.tif", wopt = list(datatype = "INT1U"))

I got this far then kind of got lost/confused and thought I would reach out for some help:

import xarray
import rioxarray

data = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]  # Example data
coords = {'x': [0, target_res, target_res*2], 'y': [0, target_res, target_res*2]} 
ds = xr.Dataset({'data': (['y', 'x'], data)}, coords=coords)
ds.rio.set_spatial_dims(x_dim='x', y_dim='y', inplace=True)
ds.rio.write_crs("EPSG:4326", inplace=True)
ds.rio.write_transform(inplace=True)
with rasterio.open(
        "test.tif, 'w',
        driver='GTiff',
        transform = ds.rio.transform(),
        crs=ds.rio.crs,
        dtype=rasterio.float32,
        #nodata=ds.rio.nodata,
        count=1, 
        width=ref.rio.width,
        
        height=ref.rio.height) as dst:
    dst.write(ds.values)

Solution

  • Part of your confusion is in mixing up your use of rasterio and rioxarray. Though rioxarray extends xarray to use rasterio functionality, it's best to either open things with rasterio or rioxarray, to perform a particular process.

    I actually do something similar in many of my workflows.

    1. Open a "mask raster" with rioxarray. You could just as easily use some other means of generating your "template" array. But make it an xarray dataset. ...

      Open and format mask raster

      print('Starting mask')
      mask_array = rioxr.open_rasterio(self.masking_raster, band_as_variable=True, masked=True, chunks=self.chunk_return)
      if not mask_array.rio.crs:
          mask_array.rio.set_crs(self.crs, inplace=True)
          mask_array.rio.write_crs(self.crs, inplace=True)
      mask_crs = mask_array.rio.crs
      print(f"Mask CRS: {mask_crs}")
      mask_bounds = mask_array.rio.bounds()
      mask_ds = mask_array.rename_vars({"band_1": "mask_values"})
      

    Since I often use data-dense input rasters as my "template" arrays/datasets, I will simplify it by stripping the actual values in the dataset variables first. ...

    # Simplify mask dataset
    if mask_ds.nbytes / 1e+9 > 90:
        mask_ds = self.mask_ds_values(mask_ds, "mask_values")
    print(f'Mask: \n----\n{mask_ds}')
    
    1. Then, open the rasters to "perform work" on, also as xarray datasets:

      Open a target/analysis raster

      grid_array = rioxr.open_rasterio(self.raw_raster, band_as_variable=True)
      grid_array.rio.write_crs(self.crs, inplace=True)
      grid_ds = grid_array.rename_vars({"band_1": "wse"})
      if not grid_ds.rio.crs:
          grid_ds.rio.set_crs(self.crs, inplace=True)
      
    2. Perform alignment using reproject.mask: ...

      grid_ds = grid_ds.rio.reproject_match(mask_ds)
      
    3. Finally, use a specific type of xarray concatenate functionality, called combine_by_coords, so that I actually hold all of my analysis arrays, as well as my mask array, as dataset variables, in a single xarray dataset. This ensures alignment and easy processing afterward:

      to_combine = to_combine + [grid_ds, mask_ds]
      comb_ds = xr.combine_by_coords(to_combine, compat="no_conflicts",
                                     combine_attrs="drop_conflicts")
      comb_ds = comb_ds.chunk(self.chunk_return)