Create a raster stack in from xarray dataset in Python

Viewed 374

I am trying to make a raster stack from an xarray dataset which I obtained from multiple netCDF files. There are 365 netCDF files each containing a 2D data of Sea Surface Temperature (SST) having height and width 3600 and 7200 respectively. To perform further operations I need to prepare a raster stack.

import netCDF4 as nc
import rasterio as rio
import numpy as np
import xarray as xr
import os

fpath = '/home/sst/2015'
pattern = '*.nc'
filelist = []


for root, dirs, files in os.walk(fpath):
    for name in files:
        filelist.append(os.path.join(fpath,name))

ds = xr.open_mfdataset(eflist, concat_dim='time', parallel = True) # netCDF data

ds_data = ds.sel(time='2015')['SST']   # xarray dataset with dimension 365x3600x7200.  

Raster stack of this xarray will be used to extract data values at point locations. I am currently using numpy and rasterio as mentioned in rasterio documentation. By iterating over the 3D xarray following code writes 365 files to HDD and later I can read and stack those.

from rasterio.transform import from_origin
transform = from_origin(-180,90, 0.05, 0.05)

fpath = '/home/sst/sst_tif'
fname = 'sst_array_'
extname = '.tiff'
timedf = ds.time    # time dimension to loop over

for i in range(len(timedf)):
    np_array = np.array(ds_data[i])
    iname = str(i)
    fwname = fpath + fname + iname + extname

    sst_tif = rio.open(fwname,
                        'w',
                       driver = 'GTiff',
                       height = 3600,
                       width = 7200,
                       dtype = np_array.dtype,
                       count = 1,
                       crs = 'EPSG:4326',
                       transform = transform)
    
    sst_tif.write(np_array, 1)
    sst_tif.close()

However this takes a very long time to process entire dataset. I also attempted converting entire xarray to a numpy array and writing all 365 layers in a single file but it freezes the Python kernel.

Is there any way I can create this raster stack in memory and do further processing without having to write it to file/s on HDD. I am trying to obtain functionality similar to that of stack function available in raster package of R.

0 Answers
Related