Last active
June 21, 2017 12:05
-
-
Save kedziorm/619141108839f154a8ded009d9cb30e8 to your computer and use it in GitHub Desktop.
trying to download data using Pydap - according to https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+PyDAP
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# -*- coding: utf-8 -*- | |
################### | |
# Author: Mateusz Kędzior | |
# Downloads sample data from NASA GSFC server | |
# just execute getSampleData() | |
################### | |
# Prerequisities: | |
# -----1--------- | |
# Create your account and approve 'NASA GESDISC DATA ARCHIVE' application as described at: | |
# http://disc.sci.gsfc.nasa.gov/registration/authorizing-gesdisc-data-access-in-earthdata_login | |
# -----2--------- | |
# Create .netrc file - I wrote small gist for that which will be executed to create .netrc file (you will be asked for password) | |
# -----3--------- | |
# Install Pydap and dependencies | |
# pip install Pydap numpy singledispatch Webob Jinja2 docopt gunicorn six mechanicalsoup | |
## or even better create virtual env for pydap: | |
## sudo apt-get install virtualenv | |
## virtualenv pydap | |
## source pydap/bin/acitvate | |
## pip install pydap | |
# -----4--------- | |
# Ensure that you have the latest pip modules installed | |
# (sample error: no cas.urs module): | |
## pip install --upgrade pip | |
## pip freeze --local | grep -v '^\-e' | cut -d = -f 1 | xargs -n1 pip install -U | |
# Samples from: https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+PyDAP | |
# are **outdated** | |
# Proper samples could be found at: http://www.pydap.org/en/latest/client.html#urs-nasa-earthdata | |
######################## | |
# Create .netrc login file needed to log in | |
import requests | |
url = "https://gist.githubusercontent.com/kedziorm/3f6634f0231ad3b743461519a47d7c9b/raw" | |
r = requests.get(url).text | |
exec(r) | |
######################## | |
# Read login and password from .netrc file | |
import netrc | |
authData = netrc.netrc().hosts['urs.earthdata.nasa.gov'] | |
myLogin = authData[0] | |
myPassword = authData[2] | |
######################## | |
import cookielib | |
import urllib2 | |
import re, uuid | |
import pydap.lib | |
from pydap.exceptions import ClientError | |
import logging | |
log = logging.getLogger("__main__")###(__name__) | |
log.setLevel(logging.DEBUG) | |
from os.path import expanduser | |
import subprocess,os.path | |
home = expanduser("~") | |
plik = os.path.join(home, str(uuid.uuid4().get_hex()) + '.tif') | |
# Set the debug level for urllib2. | |
debuglevel=1 | |
def createRasterFile(lon, lat, array, outputFile = 'myraster.tif', noData = -9.99000000e+08, myFormat='GTiff'): | |
# Saves numpy.ndarray to GDAL file | |
# Sample formats: GTiff (GeoTIFF) | |
# (AAIGrid (ASCII GRID) not supported!!) | |
import numpy as np | |
from osgeo import gdal | |
from osgeo import gdal_array | |
from osgeo import osr | |
# Remove single-dimensional entries | |
# for example change matrix of following dimenstions: (1,1,24,34) to (24,34) | |
array = np.squeeze(array) | |
# set NoData Values - I think that it is better to set in different way | |
# so commenting following line and do that by SetNoDataValue | |
# array[array==noData] = np.nan | |
xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()] | |
nrows,ncols = np.shape(array) | |
xres = (xmax-xmin)/float(ncols) | |
yres = (ymax-ymin)/float(nrows) | |
geotransform=(xmin,xres,0,ymax,0, -yres) | |
print "xmin,ymin,xmax,ymax are: %d, %d, %d, %d" % (xmin,ymin,xmax,ymax) | |
print "Saving '%s'" % outputFile | |
output_raster = gdal.GetDriverByName(myFormat).Create(outputFile,ncols, nrows, 1 ,gdal.GDT_Float32) # Open the file | |
output_raster.SetGeoTransform(geotransform) # Specify its coordinates | |
output_raster.GetRasterBand(1).SetNoDataValue(noData) | |
srs = osr.SpatialReference() # Establish its coordinate encoding | |
srs.ImportFromEPSG(4326) # This one specifies WGS84 lat long. | |
output_raster.SetProjection( srs.ExportToWkt() ) | |
output_raster.GetRasterBand(1).WriteArray(array) | |
def tif2ASCII(inFile): | |
import os | |
newFile = os.path.splitext(inFile)[0]+'.asc' | |
os.system('gdal_translate -of AAIGrid "%s" "%s"' % (inFile, newFile)) | |
def tif2XYZ(inFile): | |
# TODO: Something wrong with NA (noData) value | |
import os | |
newFile = os.path.splitext(inFile)[0]+'.txt' | |
os.system('gdal2xyz.py "%s" "%s"' % (inFile, newFile)) | |
def getSampleData(): | |
global myLogin | |
global myPassword | |
import urllib,sys | |
modules = ("Pydap", "numpy", "singledispatch", "Webob", "Jinja2", "docopt", "gunicorn", "six", "mechanicalsoup") | |
if not ('pydap' in sys.modules): | |
print("Please install Pydap and all required modules:") | |
print(modules) | |
return | |
# In case of previous query, I received and 'Exception: Unable to parse token: <!DOCTYPE' | |
# and dds had following content: | |
# | |
# <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"> | |
# <html><head> | |
# <title>301 Moved Permanently</title> | |
# </head><body> | |
# <h1>Moved Permanently</h1> | |
# <p>The document has moved <a href="https://hydro1.gesdisc.eosdis.nasa.gov/dods/_expr_{GLDAS_NOAH025SUBP_3H}{ave(rainf,time=00Z29Sep2016,time=00Z30Sep2016)}{17.00:25.25,48.75:54.50,1:1,00Z29Sep2016:00Z29Sep2016}.ascii.dds?result">here</a>.</p> | |
# </body></html> | |
# | |
# old_query = 'http://hydro1.sci.gsfc.nasa.gov/dods/_expr_%7BGLDAS_NOAH025SUBP_3H%7D%7Bave(rainf,time=00Z29Sep2016,time=00Z30Sep2016)%7D%7B17.00:25.25,48.75:54.50,1:1,00Z29Sep2016:00Z29Sep2016%7D.ascii?result' | |
# | |
sample_query = 'https://hydro1.gesdisc.eosdis.nasa.gov/dods/_expr_{GLDAS_NOAH025SUBP_3H}{ave(rainf,time=00Z29Sep2016,time=00Z30Sep2016)}{17.00:25.25,48.75:54.50,1:1,00Z29Sep2016:00Z29Sep2016}.ascii.dds?result' | |
# | |
# As a result of changing query, now I'm receiving different error: | |
# ValueError: need more than 1 value to unpack | |
# | |
from pydap.client import open_url | |
from pydap.cas.urs import setup_session | |
session = setup_session(myLogin, myPassword) | |
dataset = open_url(sample_query, session=session) | |
# Following line prints on the screen: | |
# <BaseType with data BaseProxy('https://hydro1.gesdisc.eosdis.nasa.gov/dods/_expr_{GLDAS_NOAH025SUBP_3H}{ave(rainf,time=00Z29Sep2016,time=00Z30Sep2016)}{17.00:25.25,48.75:54.50,1:1,00Z29Sep2016:00Z29Sep2016}.ascii.dds', 'result.result', dtype('>f4'), (1, 1, 24, 34), (slice(0, 1, 1), slice(0, 1, 1), slice(0, 24, 1), slice(0, 34, 1)))> | |
print(dataset.result.result) | |
# Following line causing ValueError: need more than 1 value to unpack | |
array = dataset.result.result[:] | |
# Line above cause: | |
# ValueError Traceback (most recent call last) | |
# <ipython-input-12-fca8c5905553> in <module>() | |
# ----> 1 getSampleData() | |
# <ipython-input-11-ff09df6f7236> in getSampleData() | |
# 144 dataset = open_url(sample_query, session=session) | |
# 145 print(dataset.result.result) | |
# --> 146 array = dataset.result.result[:] | |
# 147 array_dim = array.shape | |
# 148 lat = dataset.result.lat[:] | |
# /usr/local/lib/python2.7/dist-packages/pydap/model.pyc in __getitem__(self, index) | |
# 256 return np.vectorize(decode_np_strings)(self.data[index]) | |
# 257 else: | |
# --> 258 return self.data[index] | |
# 259 | |
# 260 def __len__(self): | |
# /usr/local/lib/python2.7/dist-packages/pydap/handlers/dap.pyc in __getitem__(self, index) | |
# 129 r = GET(url, self.application, self.session) | |
# 130 raise_for_status(r) | |
# --> 131 dds, data = r.body.split(b'\nData:\n', 1) | |
# 132 dds = dds.decode(r.content_encoding or 'ascii') | |
# 133 | |
# ValueError: need more than 1 value to unpack | |
array_dim = array.shape | |
lat = dataset.result.lat[:] | |
lon = dataset.result.lon[:] | |
createRasterFile(lon,lat,array,plik) | |
tif2ASCII(plik) | |
tif2XYZ(plik) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hello @kedziorm!
Did you work manage to work a way around this?
I am trying to get some NASA data via pyDAP and getting the exact same error.
Thanks,
Tomas