Get data from eartH2Observe OpeNDAP server

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:

ipython notebook --pylab=inline

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.

Catalog

http://wci.earth2observe.eu/thredds/catalog.html



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')

Extract data for a catchment.



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 [ ]: