Loop through all timesteps of a netcdf in python, store calculation output in array and merge all ar-CodePudding

I have a netcdf file containing u and v components of wind. The netcdf file contains 7305 daily values. For each day I want to do some calculations, store the output in an array and then merge all arrays, so that each array will correspond to a specific day. The initial netcdf file has lats(101) and lons(129). So the final array (vort) will have (7305,101,129) dimensions. Below I have included the code I developed so far. It works if i store each day seperately in vort but does not work when I do "append". Any help?

from matplotlib import pyplot
from matplotlib.cm import get_cmap
from __future__ import print_function
from netCDF4 import Dataset,num2date,date2num
from matplotlib.colors import from_levels_and_colors
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from metpy.units import units
from datetime import datetime
from metpy.plots import StationPlot
#
import metpy.calc as mpcalc
import xarray as xr
import cartopy.crs as ccrs
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import datetime
import cartopy.feature as cfeature
import mygrads as mg
#
#
root_dir = '/users/pr007/mkaryp/'
nc = Dataset(root_dir 'ERA5_ensmean_daymean.nc')
#
lon=nc.variables['longitude'][:]
lat=nc.variables['latitude'][:]
#
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
#
u = []
v = []
vort = np.array([])
#
for i in range(7305):
    v = np.array(nc.variables['v'][i,:,:])                    
    u = np.array(nc.variables['u'][i,:,:]) 
    #   
    v = units.Quantity(v, 'm/s')
    u = units.Quantity(u, 'm/s')
    #
    # Compute dx and dy spacing for use in divergence calculation
    vort1 = np.array(mpcalc.vorticity(u, v, dx=dx, dy=dy))
    vort = np.append([vort, vort1])
    #
    print(i)

EDIT A subsample of the data file used can be found here. This netcdf file has 10 timesteps, instead of 7305.

CodePudding user response：

If I am understanding your question correctly, what you want to do is this:

lon=nc.variables['longitude'][:]
lat=nc.variables['latitude'][:]
#
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
#
u = []
v = []
vort = []
#
for i in range(10):
    v.append(np.array(nc.variables['v'][i,:,:]))                    
    u.append(np.array(nc.variables['u'][i,:,:]))
    v[i] = units.Quantity(v[i], 'm/s')
    u[i] = units.Quantity(u[i], 'm/s')
    vort.append(np.array(mpcalc.vorticity(u[i], v[i], dx=dx, dy=dy)))

Your current for loop is overwriting the results for every index i. So you didn't really store data from all days.

You can check my results using:

print(len(vort))
print(vort[0])

The length is 10, so it means it collects data from these 10 days, instead of 1 single day in your codes.

Note that my vort is a list, not a np.array, but each element inside the list is a np.array.