creating vector using netcdf into array

Question

I'm fairly new to python and have found stack overflow one of the best resources out there, now I'm hoping someone can help me with what I believe is a fairly basic question.

I'm looking to create a land mask from a list of lats and lons and rainfall data extracted from a netCDF file. I need to get the data from the netcdf file to line up so I can remove rows which have a rainfall value of '-9999.' (indicating no data because its over the ocean). I can access the file, I can create a mesh grid, but when it comes to inserting the rainfall data for the final check I'm getting odd shapes and no luck with the logical test. Can someone have a look at this code and let me know what you think?

from netCDF4 import Dataset
import numpy as np

f=Dataset('/Testing/Ensemble_grid/1970_2012_eMAST_ANUClimate_mon_evap_v1m0_197001.nc')

lat = f.variables['latitude'][:]
lon = f.variables['longitude'][:]
rainfall = np.array(f.variables['lwe_thickness_of_precipitation_amount'])
lons, lats = np.meshgrid(lon,lat)
full_ary = np.array((lats,lons))
full_lats_lons = np.swapaxes(full_ary,0,2)
rain_data = np.squeeze(rainfall,axis=(0,))
grid = np.array((full_lats_lons,rain_data))
full_grid = np.expand_dims(grid,axis=1)
full_grid_col = np.swapaxes(full_grid,0,1)
land_grid = np.logical_not(full_grid_col[:,1]==-9999.)

N1B4 · Accepted Answer · 2015-03-26 00:23:49Z

1

Here is an alternative method that simply creates a new 2D variable, landmask, where each grid cell is either 0 (ocean) or 1 (land). (I like to use 1 and 0 landmasks because you can transform it into a boolean numpy array and do quick land-averages this way.)

import netCDF4
import numpy as np

ncfile = netCDF4.('/path/to/your/ncfile.nc', 'r')
lat = ncfile.variables['lat'][:]
lon = ncfile.variables['lon'][:]
# Presuming here that rainfall is 2D, if not, just read in the first time step, i.e. [0,:,:]
rain = ncfile.variables['lwe_thickness_of_precipitation_amount'][:,:] 
ncfile.close()

nlat, nlon = len(lat), len(lon)
# Populate a 2D landmask array, where 1=land and 0=ocean
landmask = np.zeros([nlat, nlon], dtype='int')
for y in range(nlat):
    for x in range(nlon):
        if rain[y,x]!=-9999: # We're at a land point
             landmask[y,x] = 1

# Now you can write out the landmask into a new netCDF file
filename_out = './landmask.nc'
ncfile_out = netCDF4.Dataset(filename_out, 'w')
ncfile_out.createDimension('lat', nlat)
ncfile_out.createDimension('lon', nlon)
lat_out = ncfile_out.createVariable('lat', 'f4', ('lat',))
lon_out = ncfile_out.createVariable('lon', 'f4', ('lon',))
landmask_out = ncfile_out.createVariable('landmask', 'i', ('lat', 'lon',))
setattr(lat_out, 'units', 'degrees_north')
setattr(lon_out, 'units', 'degrees_east')
setattr(landmask_out, 'description', '1=land 0=ocean')
lat_out[:] = lat
lon_out[:] = lon
landmask_out[:,:] = landmask[:,:]
ncfile_out.close()

answered Mar 26, 2015 at 0:23

N1B4

3,4631 gold badge23 silver badges25 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

IanM_63 Over a year ago

OK, thanks for the detailed response. This is substantially more involved than I expected. I was just looking to create the 2D array with only land rows (i.e.: where rainfall !=-9999.), however, this might be more suitable for the next step where I have to use the landmask to extract files. Thanks again!

Mike Laren · Accepted Answer · 2015-03-26 00:25:23Z

0

Ian, you need to put a repeatable example up here...

I suspect what you need is something like this;

x = np.array([[1, 2, 3], [4, 5, 6]], np.int32)
x.flat

edited Mar 26, 2015 at 0:25

Mike Laren

8,23618 gold badges53 silver badges72 bronze badges

answered Mar 26, 2015 at 0:22

belgaroo

1

1 Comment

IanM_63 Over a year ago

Yes, thanks Mike, I totally agree. This is my first question on stack overflow, I'll remember this next time. As to your suggestion, this does create the vector, but only creates an iterative which seemed to cause problems when I then had to search and remove rows in the array with rainfall = -9999. I ended using .flatten and that seemed to do the trick. Thanks again!

Collectives™ on Stack Overflow

creating vector using netcdf into array

2 Answers 2

1 Comment

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related