In [2]:
## Tutorial on how to work with different flavours of web services for Earth Obersvation data
# 1. netCDF-CF-OPeNDAP using netcdf4 package
# 2. OGC WCS using owslib package
# 3. OGC WMS using owslib package
# For overview see OpenEarrth.eu > Data
from IPython.display import Image
Image(url='https://publicwiki.deltares.nl/download/attachments/42401809/OpenEarthBuildingBlocks_standards_client_server.png')
Out[2]:
In [3]:
from netCDF4 import Dataset # includes OPeNDAP library
In [4]:
# 1. OPeNDAP
# define url and open dataset
url = 'http://thredds.jpl.nasa.gov/thredds/dodsC/ncml_aggregation/Chlorophyll/seawifs/aggregate__SEAWIFS_L3_CHLA_MONTHLY_9KM_R.ncml'
ds = Dataset(url)
In [5]:
# get time
print(ds.variables['time'][0],ds.variables['time'].units)
from netCDF4 import num2date
time = num2date(ds.variables['time'][:], units=ds.variables['time'].units)
print(time[0],time[1])
In [6]:
# get full coordinate sticks so be able to make geographical selection
lon = ds.variables['lon'][:]
lat = ds.variables['lat'][:]
[lon[0],lat[0],lon[-1],lat[-1]]
Out[6]:
In [7]:
# from bbox get netCDF-OPeNDAP subsetting indices
import numpy as np
bbox = (-10., 50., 15., 60.)
ilon = []
for i,il in enumerate(lon):
if (il >= bbox[0]) & (il <= bbox[2]):
ilon.append(i)
ilat = []
for i,il in enumerate(lat):
if (il >= bbox[1]) & (il <= bbox[3]):
ilat.append(i)
[ilon[0],ilat[0],ilon[-1],ilat[-1]]
Out[7]:
In [8]:
it = 0 # first
l3m_data = ds.variables['l3m_data'][it,ilat[0]:ilat[-1],ilon[0]:ilon[-1]]
In [9]:
# data is numpy masked array
l3m_data
Out[9]:
In [10]:
clim = [0, 10]
print(clim)
[np.nanmin(l3m_data),np.max(l3m_data)]
Out[10]:
In [11]:
# close dataset
ds.close()
In [12]:
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
fig=plt.figure()
fig.set_figwidth(10)
ax = plt.axes([0.1,.8,0.8,0.8])
dd = 8
h = plt.pcolor(lon[ilon[0]:ilon[-1]],lat[ilat[0]:ilat[-1]],l3m_data)
h.set_clim([0, 2])
plt.colorbar()
plt.title(str(time[it]))
Out[12]:
In [13]:
# 2. WCS
#see also: http://stommel.tamu.edu/~baum/wcs.html
import owslib.wcs
url = 'http://thredds.jpl.nasa.gov/thredds/wcs/ncml_aggregation/Chlorophyll/seawifs/aggregate__SEAWIFS_L3_CHLA_MONTHLY_9KM_R.ncml?'
wcs = owslib.wcs.WebCoverageService(url)
In [14]:
wcs.contents
Out[14]:
In [15]:
meta = wcs.contents['l3m_data']
print(meta.boundingBoxWGS84)
print(meta.grid)
print(meta.supportedFormats)
In [16]:
bbox = (-10., 50., 15., 60.)
In [17]:
help(wcs.getCoverage)
output = wcs.getCoverage(identifier='l3m_data',
bbox=bbox,
format='GeoTIFF',
crs='OGC:CRS84',
version="1.0.0")
# save to file
f=open('wcs_test.tiff','wb')
f.write(output.read())
f.close()
In [18]:
# no PIL for python 3.4.*
from PIL import Image
im = Image.open('wcs_test.tiff')
%matplotlib inline
plt.imshow(img,extent=[bbox[0],bbox[1],bbox[2],bbox[3]])
In [19]:
# 3. WMS
# see also: http://nbviewer.ipython.org/github/Unidata/tds-python-workshop/blob/master/wms_sample.ipynb
import owslib.wms
url = 'http://thredds.jpl.nasa.gov/thredds/wms/ncml_aggregation/Chlorophyll/seawifs/aggregate__SEAWIFS_L3_CHLA_MONTHLY_9KM_R.ncml?'
wms = owslib.wms.WebMapService(url)
In [20]:
# let's see what the server offers.
wms.contents
wms.contents['l3m_data'].timepositions[0], wms.contents['l3m_data'].timepositions[-1]
Out[20]:
In [21]:
dir(wms.contents['l3m_data'])
Out[21]:
In [22]:
bbox = (-10., 50., 15., 60.)
f = wms.getmap(['l3m_data'], \
srs='EPSG:4326', \
bbox=bbox, \
size=(1024, 1024), \
format='image/png', \
time='1997-09-01T00:00:00.000Z',colorscalerange='0,10')
# colorscalerange does not work
In [23]:
# convert binary data to IO stream
import matplotlib.pyplot as plt
from io import BytesIO, StringIO
f_io = BytesIO(f.read())
img = plt.imread(f_io)
# x = bbox[0] + (bbox[2]-bbox[0])*np.arange(0,img.shape[0]+1)/img.shape[0]
# y = bbox[1] + (bbox[3]-bbox[1])*np.arange(0,img.shape[1]+1)/img.shape[1]
In [24]:
%matplotlib inline
plt.imshow(img,extent=[bbox[0],bbox[1],bbox[2],bbox[3]])
Out[24]:
In [ ]: