Estimate dataset size from the OPeNDAP DDS. Here we use regular expressions to parse the DDS and just the variable size (32 or 64 bit Int or Float) by their shapes. This represents the size in memory, not on disk, since the data could be compressed. But the data in memory is in some sense a more true representation of the quantity of data available by the service.
In [85]:
from owslib.csw import CatalogueServiceWeb
from owslib import fes
import pandas as pd
import datetime as dt
import requests
import re
import time
In [2]:
def service_urls(records,service_string='urn:x-esri:specification:ServiceType:odp:url'):
"""
Get all URLs matching a specific ServiceType
Unfortunately these seem to differ between different CSW-ISO services.
For example, OpenDAP is specified:
NODC geoportal: 'urn:x-esri:specification:ServiceType:OPeNDAP'
NGDC geoportal: 'urn:x-esri:specification:ServiceType:odp:url'
"""
urls=[]
for key,rec in records.iteritems():
#create a generator object, and iterate through it until the match is found
#if not found, gets the default value (here "none")
url = next((d['url'] for d in rec.references if d['scheme'] == service_string), None)
if url is not None:
urls.append(url)
return urls
In [5]:
endpoint = 'http://www.ngdc.noaa.gov/geoportal/csw' # NGDC/IOOS Geoportal
csw = CatalogueServiceWeb(endpoint,timeout=60)
csw.version
Out[5]:
In [6]:
[op.name for op in csw.operations]
Out[6]:
In [12]:
for oper in csw.operations:
if oper.name == 'GetRecords':
print oper.constraints
Since the supported ISO queryables contain apiso:ServiceType, we can use CSW to find all datasets with services that contain the string "dap"
In [4]:
val = 'dap'
service_type = fes.PropertyIsLike(propertyname='apiso:ServiceType',literal=('*%s*' % val),
escapeChar='\\',wildCard='*',singleChar='?')
filter_list = [ service_type]
In [13]:
csw.getrecords2(constraints=filter_list,maxrecords=10000,esn='full')
len(csw.records.keys())
Out[13]:
By printing out the references from a random record, we see that for this CSW the DAP URL is identified by
urn:x-esri:specification:ServiceType:odp:url
In [14]:
choice=random.choice(list(csw.records.keys()))
print choice
csw.records[choice].references
Out[14]:
Get all the DAP endpoints
In [15]:
dap_urls = service_urls(csw.records,service_string='urn:x-esri:specification:ServiceType:odp:url')
len(dap_urls)
Out[15]:
In [17]:
def calc_dsize(txt):
'''
Calculate dataset size from the OPeNDAP DDS.
Approx method: Multiply 32|64 bit Int|Float variables by their shape.
'''
# split the OpenDAP DDS on ';' characters
all = re.split(';',txt)
'''
Use regex to find numbers following Float or Int (e.g. Float32, Int64)
and also numbers immediately preceding a "]". The idea is that in line like:
Float32 Total_precipitation_surface_6_Hour_Accumulation[time2 = 74][y = 303][x = 491];
we want to find only the numbers that are not part of a variable or dimension name
(want to return [32, 74, 303, 491], *not* [32, 6, 2, 74, 303, 491])
'''
m = re.compile('\d+(?=])|(?<=Float)\d+|(?<=Int)\d+')
dsize=0
for var in all:
c = map(int,m.findall(var))
if len(c)>=2:
vsize = reduce(lambda x,y: x*y,c)
dsize += vsize
return dsize/1.0e6/8. # return megabytes
In [25]:
def tot_dsize(url,timeout=10):
das = url + '.dds'
tot = 0
try:
response = requests.get(das,verify=True, timeout=timeout)
except:
return tot, -1
if response.status_code==200:
# calculate the total size for all variables:
tot = calc_dsize(response.text)
# calculate the size for MAPS variables and subtract from the total:
maps = re.compile('MAPS:(.*?)}',re.MULTILINE | re.DOTALL)
map_text = ''.join(maps.findall(response.text))
if map_text:
map_tot = calc_dsize(map_text)
tot -= map_tot
return tot,response.status_code
In [ ]:
from __future__ import print_function
time0 = time.time()
good_data=[]
bad_data=[]
count=0
for url in dap_urls:
count += 1
dtot, status_code = tot_dsize(url,timeout=2)
if status_code==200:
good_data.append([url,dtot])
print('[{}]Good:{},{}'.format(count,url,dtot), end='\r')
else:
bad_data.append([url,status_code])
print('[{}]Fail:{},{}'.format(count,url,status_code), end='\r')
print('Elapsed time={} minutes'.format((time.time()-time0)/60.))
In [34]:
print('Elapsed time={} minutes'.format((time.time()-time0)/60.))
In [22]:
len(good_data)
Out[22]:
In [23]:
len(bad_data)
Out[23]:
In [35]:
bad_data[0][0]
Out[35]:
Loop through the datasets that failed in the 2 second timeout to see if any of them work with a 10 second timeout
In [ ]:
time0 = time.time()
good_data2=[]
bad_data2=[]
count=0
for item in bad_data:
url = item[0]
count += 1
dtot, status_code = tot_dsize(url,timeout=10)
if status_code==200:
good_data2.append([url,dtot])
print('[{}]Good:{},{}'.format(count,url,dtot), end='\r')
else:
bad_data2.append([url,status_code])
print('[{}]Fail:{},{}'.format(count,url,status_code), end='\r')
print('Elapsed time={} minutes'.format((time.time()-time0)/60.))
yipes, that took forever with 10 second timeout. How many more datasets did we get?
In [87]:
len(bad_data)-len(bad_data2)
Out[87]:
So how much data are we serving?
In [52]:
sum=0
for ds in good_data:
sum +=ds[1]
print('{} terabytes'.format(sum/1.e6))
How much more data do we get if we allow 10 second timeout instead of 2?
In [53]:
sum=0
for ds in good_data2:
sum +=ds[1]
print('{} terabytes'.format(sum/1.e6))
In [59]:
url=[]
size=[]
for item in good_data:
url.append(item[0])
size.append(item[1])
In [55]:
d={}
d['url']=url
d['size']=size
In [88]:
good = pd.DataFrame(d)
In [89]:
good.head()
Out[89]:
In [71]:
df2=df.sort(['size'],ascending=0)
In [76]:
df2.head()
Out[76]:
In [91]:
url=[]
code=[]
for item in bad_data:
url.append(item[0])
code.append(item[1])
In [92]:
d={}
d['url']=url
d['code']=code
bad = pd.DataFrame(d)
In [93]:
bad.head()
Out[93]:
In [102]:
bad.to_csv('bad.csv')
In [103]:
good.to_csv('good.csv')
In [ ]: