CODEWITHSUNDEEP
pythonnetcdfnetcdf4
Muon
Muon

Reputation: 1346

Python - Use list of points to extract data from gridded NetCDF without for loops

The following example uses the "Unidata" sample netCDF dataset of eastward wind which can be downloaded from here (2.8 MB)

I have two lists of integers that correspond to the x and y index of a gridded array in a netCDF file. I want to extract the data and save it to a 1 dimensional array or list for each of the point combinations (e.g. points: [(x[0],y[0]), (x[1],y[1]), (x[2],y[2]), ... , (x[n],y[n])]).

I can do this quite easily using this method...

from netCDF4 import Dataset

# grid point lists
lat = [20, 45, 56, 67, 88, 98, 115]
lon = [32, 38, 48, 58, 87, 92, 143]

# open netCDF file
nc_file = "./sresa1b_ncar_ccsm3-example.nc"
fh = Dataset(nc_file, mode='r')

# extract variable
point_list = zip(lat,lon)
ua_list = []
for i, j in point_list:
    ua_list.append(fh.variables['ua'][0,16,i,j])

print(ua_list)

Which returns:

[59.29171, 17.413916, -4.4006901, -11.15424, -5.2684789, 2.1235929, -6.134573]

However append() is clunky on big datasets and I'm trying to speed up my code so I also do not want to use a for loop and would rather return the results in a single line. I've tried doing so using this line:

# extract variable
ua_array = fh.variables['ua'][0,16,lat,lon]
print(ua_array)

Which returns every single possible combination of points instead of just the ones I am after:

[[ 59.2917099   60.3418541   61.81352234  62.66215515  60.6419754 60.00745392  52.48550797]
[ 18.80122566  17.41391563  14.83201313  12.67425823  13.99616718 14.4371767   14.12419605]
[ -5.56457043  -5.20643377  -4.40069008  -3.25902319  -2.36573601 -2.25667071  -1.0884304 ]
[-11.66207981 -11.46785831 -11.35252953 -11.15423965 -11.35271263 -11.55139542 -11.68573093]
[ -1.15064895  -1.52471519  -2.12152767  -2.67548943  -5.26847887 -5.79328251  -6.16713762]
[ -1.95770085  -0.56232995   0.82722098   1.39629912   2.65125418 2.12359285  -6.47501516]
[ -9.76508904 -10.13490105 -10.76805496 -11.31607246 -11.93865585 -11.56440639  -6.13457298]]

How can I slice the netCDF file so I can get the same result as the above code in a single line? Thanks in advance.

Upvotes: 3

Views: 2702

Answers (1)

Mike Müller
Mike Müller

Reputation: 85462

Do normal indexing with 0and 16 first, followed by advanced indexing with lat and lon: 

ua_array = fh.variables['ua'][0,16][lat,lon]
print(ua_array)

Output:

[ 59.2917099   17.41391563  -4.40069008 -11.15423965  -5.26847887
   2.12359285  -6.13457298]

BTW, ua_array is a NumPy array. Therefore, calling its ua_list is bit misleading.

Upvotes: 3

Related Questions