0

I'm preprocessing HMASR data, the data contains multiple data files with format of netcdf, and each tile (1X1 degree) is a nc file. So I read it through xarray open_mfdataset() to mosaic it spatially and do transformation after this. The data after mosaicing is like 4050x9900x366 which is listed below:

<xarray.Dataset>
Dimensions:    (Latitude: 4050, Longitude: 9900, Day: 366, Stats: 5)
Coordinates:
  * Latitude   (Latitude) float64 27.0 27.01 27.01 27.02 ... 44.99 44.99 45.0
  * Longitude  (Longitude) float64 61.0 61.01 61.01 61.02 ... 105.0 105.0 105.0
Dimensions without coordinates: Day, Stats
Data variables:
    SWE_Post   (Day, Stats, Longitude, Latitude) float32 dask.array<chunksize=(366, 5, 225, 225), meta=np.ndarray>
    SCA_Post   (Day, Stats, Longitude, Latitude) float32 dask.array<chunksize=(366, 5, 225, 225), meta=np.ndarray>

The speed is acceptable so far, but it becomes desperately slow when it comes to write it to the geotif.

The raster I want to write shaped (Latitude: 4050, Longitude: 9900). Firstly, I tried to use xarray function to_raster to write it , and with resample to match the template(you can see it from the commented code), I found that the writing processes itself are slow enough, so I gave it up. It kept getting performance warning. I also tried the convert the data to numpy.array to output it with gdal, but it gives no improvement. I think the problem is linked with the dask or chunksize. As i have seen example of the larger size image writing without performance issue.

This is part of my code:

import numpy as np
import xarray as xr
from lytools import * # a customed package
from tqdm import tqdm
from osgeo import gdal
import rioxarray
import dask    
def SWE_merged_annually(self, indir):
    #### Below is code used to create the array that will be written to a GeoTIFF ####
    # mosaic
    ds = xr.open_mfdataset(indir,
                           combine="by_coords")
    SWE_year = ds.SWE_Post.sel(Stats=0)
    SWE_peak_2000 = xr.apply_ufunc(lambda x: x.max(axis=-1),
                                   SWE_year,
                                   input_core_dims=[['Day']],
                                   dask='allowed',
                                   vectorize=True)

    SWE_peak_dowy_2000 = xr.apply_ufunc(lambda x: x.argmax(axis=-1),
                                        SWE_year,
                                        input_core_dims=[['Day']],
                                        dask='allowed',
                                        vectorize=True)
    # mask
    SWE_peak_2000 = SWE_peak_2000.where(self.classifi == 1)
    SWE_peak_2000 = SWE_peak_2000.where(self.Non_seasonal_snow == 0)
    SWE_peak_dowy_2000 = SWE_peak_dowy_2000.where(self.classifi == 1)
    SWE_peak_dowy_2000 = SWE_peak_dowy_2000.where(self.Non_seasonal_snow == 0)

    ##### Above is code used to create array to write to GeoTIFF #####
    ##### writing part ######
    outdir1 = "/Volumes/My Passport/组内数据/雪水当量/HMA_SR/SWE_peak_tiff/"
    outdir2 = "/Volumes/My Passport/组内数据/雪水当量/HMA_SR/SWE_peal_DOWY_tiff/"
    Tools().mkdir(outdir1)
    Tools().mkdir(outdir2)

    outfname1 = outdir1 + str(self.year_now) + '.tif'
    outfname2 = outdir2 + str(self.year_now) + '.tif'

    # scheme 1 writing with rioxarray
    # ds = rioxarray.open_rasterio(tif_template_yuan, masked=True).isel(band=0)
    # SWE_peak_2000_resample = SWE_peak_2000.squeeze().transpose('Latitude', 'Longitude').rio. \
    #         set_spatial_dims(x_dim='Longitude', y_dim='Latitude').rio.write_crs('epsg:4326').rio.reproject_match(ds,resampling=5) 
    # SWE_peak_2000_dowy_resample = SWE_peak_dowy_2000.squeeze().transpose('Latitude', 'Longitude').rio. \
    #         set_spatial_dims(x_dim='Longitude', y_dim='Latitude').rio.write_crs('epsg:4326').rio.reproject_match(ds, resampling=5)
    # SWE_peak_2000_resample.to_raster(outfname1)
    # SWE_peak_2000_dowy_resample.to_raster(outfname2)

    # scheme 2 write with gdal
    Lon = SWE_peak_2000.Longitude
    Lat = SWE_peak_2000.Latitude
    LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
    SWE_peak_arr = np.transpose(np.array(SWE_peak_2000))
    SWE_peak_dowy_arr = np.transpose(np.array(SWE_peak_dowy_2000))
    ToRaster().array2raster(outfname1, LonMin, LatMax, Lon_Res, -Lat_Res,
                            SWE_peak_arr)  # this is a customed function supposed to write tiff
    ToRaster().array2raster(outfname2, LonMin, LatMax, Lon_Res, -Lat_Res, SWE_peak_dowy_arr)

The performance warnings are below, i get rid of it by add dask.config.set(**{'array.slicing.split_large_chunks': False}), but i'm not sure if it helps the progress substantially

PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

Anyone know how to improve the performance of this code. And hopefully with resample functionality.

2
  • I just finished writing 4 Geotiffs to my disk, it costs me 9 hours. I think it's abnormally slow, as each Geotiff is no more than 200MB. I have similar question when dealing with CMIP6 data, it's also about xarray and dask, I think it's conneted. I tried add some configuration to dask, it's also won't help.
    – 孟泽楷
    Commented Mar 5 at 2:49
  • Nevermind, i finished runing it, it's very slow though. i check some github issues, lots of question tagged with performance have quite a few comments. I guess the it's just what it is.
    – 孟泽楷
    Commented Mar 14 at 6:37

0

Your Answer

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

Browse other questions tagged or ask your own question.