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.