Solución:
Te recomiendo encarecidamente que eches un vistazo a la xarray y dask proyectos. El uso de estas poderosas herramientas le permitirá dividir fácilmente el cálculo en partes. Esto ofrece dos ventajas: puede calcular datos que no caben en la memoria y puede utilizar todos los núcleos de su máquina para obtener un mejor rendimiento. Puede optimizar el rendimiento eligiendo adecuadamente el tamaño del fragmento (consulte la documentación).
Puede cargar sus datos desde netCDF haciendo algo tan simple como
import xarray as xr
ds = xr.open_dataset(path_file)
Si desea fragmentar sus datos en años a lo largo de la dimensión de tiempo, especifique el chunks parámetro (asumiendo que la coordenada del año se llama ‘año’):
ds = xr.open_dataset(path_file, chunks={'year': 10})
Dado que las otras coordenadas no aparecen en el chunks dict, se utilizará un solo fragmento para ellos. (Vea más detalles en la documentación aquí). Esto será útil para su primer requisito, donde desea multiplicar cada año por una matriz 2D. Simplemente harías:
ds['new_var'] = ds['var_name'] * arr_2d
Ahora, xarray y dask están calculando su resultado perezosamente. Para activar el cálculo real, simplemente puede preguntar xarray para guardar su resultado en netCDF:
ds.to_netcdf(new_file)
El cálculo se activa a través de dask, que se encarga de dividir el procesamiento en trozos y, por lo tanto, permite trabajar con datos que no caben en la memoria. Además, dask se encargará de utilizar todos los núcleos de su procesador para procesar fragmentos.
los xarray y dask los proyectos aún no manejan bien situaciones en las que los fragmentos no se “alinean” bien para el cálculo paralelo. Dado que en este caso solo dividimos en la dimensión ‘año’, esperamos no tener problemas.
Si desea agregar dos archivos netCDF diferentes juntos, es tan simple como:
ds1 = xr.open_dataset(path_file1, chunks={'year': 10})
ds2 = xr.open_dataset(path_file2, chunks={'year': 10})
(ds1 + ds2).to_netcdf(new_file)
Proporcioné un ejemplo completamente funcional utilizando un conjunto de datos disponible en línea.
In [1]:
import xarray as xr
import numpy as np
# Load sample data and strip out most of it:
ds = xr.open_dataset('ECMWF_ERA-40_subset.nc', chunks = {'time': 4})
ds.attrs = {}
ds = ds[['latitude', 'longitude', 'time', 'tcw']]
ds
Out[1]:
<xarray.Dataset>
Dimensions: (latitude: 73, longitude: 144, time: 62)
Coordinates:
* latitude (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
* longitude (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
* time (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
tcw (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...
In [2]:
arr2d = np.ones((73, 144)) * 3.
arr2d.shape
Out[2]:
(73, 144)
In [3]:
myds = ds
myds['new_var'] = ds['tcw'] * arr2d
In [4]:
myds
Out[4]:
<xarray.Dataset>
Dimensions: (latitude: 73, longitude: 144, time: 62)
Coordinates:
* latitude (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
* longitude (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
* time (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
tcw (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...
new_var (time, latitude, longitude) float64 30.46 30.46 30.46 30.46 ...
In [5]:
myds.to_netcdf('myds.nc')
xr.open_dataset('myds.nc')
Out[5]:
<xarray.Dataset>
Dimensions: (latitude: 73, longitude: 144, time: 62)
Coordinates:
* latitude (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
* longitude (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
* time (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
tcw (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...
new_var (time, latitude, longitude) float64 30.46 30.46 30.46 30.46 ...
In [6]:
(myds + myds).to_netcdf('myds2.nc')
xr.open_dataset('myds2.nc')
Out[6]:
<xarray.Dataset>
Dimensions: (latitude: 73, longitude: 144, time: 62)
Coordinates:
* time (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
* latitude (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
* longitude (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
Data variables:
tcw (time, latitude, longitude) float64 20.31 20.31 20.31 20.31 ...
new_var (time, latitude, longitude) float64 60.92 60.92 60.92 60.92 ...
Compruebe la fragmentación del archivo. ncdump -s <infile> dará la respuesta. Si el tamaño del fragmento en la dimensión de tiempo es mayor que uno, debe leer la misma cantidad de años a la vez; de lo contrario, está leyendo varios años a la vez del disco y usando solo uno a la vez. ¿Qué tan lento es lento? Un máximo de pocos segundos por paso de tiempo suena razonable para una matriz de este tamaño. Dar más información sobre lo que hace con los datos más adelante puede brindarnos más orientación sobre dónde puede estar el problema.