This data recipe was tested under: Python 2.7.13 :: Anaconda 4.3.1 (x86_64). The data are sourced from NASA's National Snow and Ice Data Center Distributed Active Archive Center (NSIDC DAAC).
Step 1
import os
import matplotlib as mpl
import matplotlib.pyplot as plt from mpl_toolkits.basemap
import Basemap
import numpy as np
import csv
import urllib2
import math
Set OPeNDAP URL. NSIDC runs a Hyrax server that enables CF output for HDF5 data product. Such a server cannot handle compound data type in HDF-EOS5 Point data product. Therefore, we copied a sample file and put it on our demo Hyrax server that disabled CF Output.
For large subset of data, this code will not work for some reason. Thus, we select 10 points using OPeNDAP's constraint expression.
Neither PyDAP nor netCDF4 python module can handle OPeNDAP's array of structure. Thus, we will read data in ASCII instead of DAP2.
url = ("https://eosdap.hdfgroup.org:8989" # Server
"/opendap/hyrax/data/NASAFILES/hdf5/" # Path to data on server
"AMSR_U2_L2_Land_B01_201207232336_D.he5" # HDF-EOS5 Point file
".ascii?" # OPeNDAP request
"/HDFEOS/POINTS/AMSR-2%20Level%202%20Land%20Data/Data/" # Group Path
"Combined%20NPD%20and%20SCA%20Output%20Fields" # HDF5 Dataset
"[3700:1:3709]") # OPeNDAP constraint - select 10 points.
response = urllib2.urlopen(url)
cr = csv.reader(response)
# The first line of output is dataset name.
i = 0
# Every 3rd row is the actual value in CSV format.
j = -1
lat = np.array([], dtype=float)
lon = np.array([], dtype=float)
data = np.array([], dtype=float)
for row in cr:
if i != 0:
j = j+1
if i != 0 and (j % 3 == 2):
# Latitude
# print row[1]
lat = np.append(lat, float(row[1]))
# Longitude
# print row[2]
lon = np.append(lon, float(row[2]))
# SoilMoistureNPD
# print row[16]
data = np.append(data, float(row[16]))
i = i+1
print lat, lon, data
[-5.95511 -6.1588 -6.35505 -6.56324 -6.75285 -6.96071 -7.14491 -7.34555
-7.53102 -7.7309 ] [ 17.9583 17.9724 17.9519 17.9644 17.9538 17.9567 17.9625 17.952
17.9661 17.9547] [ 0.13084 0.122343 0.130008 0.140189 0.148062 0.155974 0.153412
0.144981 0.147109 0.150234]
Step 2
Let's plot the data using the above lat/lon/data ASCII values retrieved from OPeNDAP server.
m = Basemap(projection='cyl', resolution='l',
llcrnrlat=-90, urcrnrlat=90,
llcrnrlon=-180, urcrnrlon=180)
m.drawcoastlines(linewidth=0.5)
m.scatter(lon, lat, c=data, s=1, cmap=plt.cm.jet, edgecolors=None, linewidth=0)
fig = plt.gcf()
plt.show()