This small notebook shows how the get data from the eartH2Observe OpeNDAP server using python and the netCDF4 package. To run it make sure you notebook server is started like this:
If you want to run this outside of ipython just download the .py file. IN that case make sure to import numpy and scipy first.
See also: https://publicwiki.deltares.nl/display/OET/OPeNDAP+subsetting+with+python
In [ ]:
#!/usr/bin/env python
# we need the netcdf4 library. This allows us to read netcdf from OpeNDAP or files
import netCDF4
# Uncomment these if needed
#from numpy import *
#from scipy import *
# specify an url
rainurl = 'http://wci.earth2observe.eu/thredds/dodsC/ecmwf/met_forcing_v0/2012/Rainf_daily_E2OBS_201201.nc'
# create a dataset object
raindataset = netCDF4.Dataset(rainurl)
%pylab inline
In [ ]:
def getlat(ncdataset):
"""
"""
for a in ncdataset.variables:
if ncdataset.variables[a].standard_name == 'latitude':
return ncdataset.variables[a]
return None
def getlon(ncdataset):
"""
"""
for a in ncdataset.variables:
if ncdataset.variables[a].standard_name == 'longitude':
return ncdataset.variables[a]
return None
def getheigth(ncdataset):
"""
"""
for a in ncdataset.variables:
if ncdataset.variables[a].standard_name == 'height':
return ncdataset.variables[a]
return None
a= getlon(raindataset)
In [ ]:
a.shape
In [ ]:
# Get some data ......
# note that with larg set you do not want to read all in memory at once like we do here
precip = raindataset.variables['Rainf'][:]
lat = raindataset.variables['lat'][:]
lon = raindataset.variables['lon'][:]
# Show the max value in the matrix
precip[1,:,:].max()
In [ ]:
# Display the units
raindataset.variables['Rainf'].units
Now let's make a plot to see the contents of the file
In [ ]:
f, ax = plt.subplots(figsize=(16, 9))
Lon,Lat = meshgrid(lon, lat)
mesh = pcolormesh(Lon,Lat,log( 0.000001 + precip[1,:,:]),cmap='hot')
In [ ]:
# Select lat and long for our cathcment
# Bounding box for our catchment
BB = dict(
lon=[ 143, 150],
lat=[-37, -33]
)
(latidx,) = logical_and(lat >= BB['lat'][0], lat < BB['lat'][1]).nonzero()
(lonidx,) = logical_and(lon >= BB['lon'][0], lon < BB['lon'][1]).nonzero()
# get rid of the non used lat/lon now
lat = lat[latidx]
lon = lon[lonidx]
# Now get the time for the x-axis
time = raindataset.variables['time']
timeObj = netCDF4.num2date(time[:], units=time.units, calendar=time.calendar)
In [ ]:
#Now determine area P for each timestep and display in a graph
# first the mean per area lat, next average those also
# Multiply with timestep in seconds to get mm
p_select = precip[:,latidx.min():latidx.max(),lonidx.min():lonidx.max()]
In [ ]:
# PLot the sum over this month for the subcatchment
Lon,Lat = meshgrid(lon, lat)
mesh = pcolormesh(Lon,Lat,p_select.sum(axis=0))
title("Cumulative precipitation")
In [ ]:
p_mean = p_select.mean(axis=2).mean(axis=1) * hstack((diff(time[:]),mean(diff(time[:]))))
In [ ]:
f, ax = plt.subplots(figsize=(16, 6))
ax.plot(timeObj,p_mean)
title("Precipitation")
In [ ]:
#Save to disk
p_select = precip[:,latidx.min():latidx.max(),lonidx.min():lonidx.max()]
In [ ]: