Use CSW to find the total amount of opendap model data served by NGDC

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

Find model results at NGDC


In [3]:
endpoint = 'http://www.ngdc.noaa.gov/geoportal/csw' #  NGDC/IOOS Geoportal
csw = CatalogueServiceWeb(endpoint,timeout=60)
csw.version


Out[3]:
'2.0.2'

In [4]:
[op.name for op in csw.operations]


Out[4]:
['GetCapabilities',
 'DescribeRecord',
 'GetRecords',
 'GetRecordById',
 'Transaction']

In [5]:
for oper in csw.operations:
    if oper.name == 'GetRecords':
        print oper.constraints


[Constraint: SupportedCommonQueryables - ['Subject', 'Title', 'Abstract', 'AnyText', 'Format', 'Identifier', 'Modified', 'Type', 'BoundingBox'], Constraint: SupportedISOQueryables - ['apiso:Subject', 'apiso:Title', 'apiso:Abstract', 'apiso:AnyText', 'apiso:Format', 'apiso:Identifier', 'apiso:Modified', 'apiso:Type', 'apiso:BoundingBox', 'apiso:CRS.Authority', 'apiso:CRS.ID', 'apiso:CRS.Version', 'apiso:RevisionDate', 'apiso:AlternateTitle', 'apiso:CreationDate', 'apiso:PublicationDate', 'apiso:OrganizationName', 'apiso:HasSecurityConstraints', 'apiso:Language', 'apiso:ResourceIdentifier', 'apiso:ParentIdentifier', 'apiso:KeywordType', 'apiso:TopicCategory', 'apiso:ResourceLanguage', 'apiso:GeographicDescriptionCode', 'apiso:Denominator', 'apiso:DistanceValue', 'apiso:DistanceUOM', 'apiso:TempExtent_begin', 'apiso:TempExtent_end', 'apiso:ServiceType', 'apiso:ServiceTypeVersion', 'apiso:Operation', 'apiso:OperatesOn', 'apiso:OperatesOnIdentifier', 'apiso:OperatesOnName', 'apiso:CouplingType'], Constraint: AdditionalQueryables - ['apiso:Degree', 'apiso:AccessConstraints', 'apiso:OtherConstraints', 'apiso:Classification', 'apiso:ConditionApplyingToAccessAndUse', 'apiso:Lineage', 'apiso:ResponsiblePartyRole', 'apiso:ResponsiblePartyName', 'apiso:SpecificationTitle', 'apiso:SpecificationDate', 'apiso:SpecificationDateType']]

Since the supported ISO queryables contain apiso:ServiceType, we can use CSW to find all datasets with services that contain the string "dap"


In [6]:
val = 'dap'
service_type = fes.PropertyIsLike(propertyname='apiso:ServiceType',literal=('*%s*' % val),
                        escapeChar='\\',wildCard='*',singleChar='?')
filter_list = [ service_type]

In [7]:
csw.getrecords2(constraints=filter_list,maxrecords=10000,esn='full')
len(csw.records.keys())


Out[7]:
2137

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 [8]:
choice=random.choice(list(csw.records.keys()))
print choice
csw.records[choice].references


buoy43176-agg
Out[8]:
[{'scheme': 'urn:x-esri:specification:ServiceType:distribution:url',
  'url': 'http://sos.maracoos.org/stable/dodsC/sldmb/buoy43176-agg.ncml.html'},
 {'scheme': 'urn:x-esri:specification:ServiceType:distribution:url',
  'url': 'http://www.ncdc.noaa.gov/oa/wct/wct-jnlp-beta.php?singlefile=http://sos.maracoos.org/stable/dodsC/sldmb/buoy43176-agg.ncml'},
 {'scheme': 'urn:x-esri:specification:ServiceType:sos:url',
  'url': 'http://sos.maracoos.org/stable/sos/sldmb/buoy43176-agg.ncml?service=SOS&version=1.0.0&request=GetCapabilities'},
 {'scheme': 'urn:x-esri:specification:ServiceType:odp:url',
  'url': 'http://sos.maracoos.org/stable/dodsC/sldmb/buoy43176-agg.ncml'},
 {'scheme': 'urn:x-esri:specification:ServiceType:download:url',
  'url': 'http://sos.maracoos.org/stable/dodsC/sldmb/buoy43176-agg.ncml.html'}]

Get all the DAP endpoints


In [9]:
dap_urls = service_urls(csw.records,service_string='urn:x-esri:specification:ServiceType:odp:url')
len(dap_urls)


Out[9]:
2027

In [10]:
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 [11]:
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 [12]:
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.))


Elapsed time=16.6359228333 minutes

In [13]:
print('Elapsed time={} minutes'.format((time.time()-time0)/60.))


Elapsed time=16.6360023499 minutes

In [14]:
len(good_data)


Out[14]:
1559

In [15]:
len(bad_data)


Out[15]:
468

In [16]:
bad_data[0][0]


Out[16]:
'http://www.neracoos.org/thredds/dodsC/UMO/DSG/SOS/A01/Doppler/HistoricRealtime/Agg.ncml'

Loop through the datasets that failed in the 2 second timeout to see if any of them work with a 10 second timeout


In [17]:
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.))


Elapsed time=34.2986972332 minutes

yipes, that took forever with 10 second timeout. How many more datasets did we get?


In [18]:
len(bad_data)-len(bad_data2)


Out[18]:
74

So how much data are we serving?


In [19]:
sum=0
for ds in good_data:
    sum +=ds[1]
    
print('{} terabytes'.format(sum/1.e6))


49.9728701633 terabytes

How much more data do we get if we allow 10 second timeout instead of 2?


In [20]:
sum=0
for ds in good_data2:
    sum +=ds[1]
    
print('{} terabytes'.format(sum/1.e6))


3.92371649923 terabytes

In [21]:
url=[]
size=[]
for item in good_data:
    url.append(item[0])
    size.append(item[1])

In [22]:
d={}
d['url']=url
d['size']=size

In [23]:
good = pd.DataFrame(d)

In [24]:
good.head()


Out[24]:
size url
0 161.589067 http://oos.soest.hawaii.edu/thredds/dodsC/paci...
1 274.258484 http://thredds.coastal.ufl.edu:8080/thredds/do...
2 721.913524 http://thredds.coastal.ufl.edu:8080/thredds/do...
3 79.214080 http://thredds.coastal.ufl.edu:8080/thredds/do...
4 208.505984 http://thredds.coastal.ufl.edu:8080/thredds/do...

In [25]:
good_sorted=good.sort(['size'],ascending=0)

In [26]:
good_sorted.head()


Out[26]:
size url
222 8075141.971392 http://ecowatch.ncddc.noaa.gov/thredds/dodsC/n...
184 7494709.539768 http://geoport.whoi.edu/thredds/dodsC/coawst_4...
325 4567598.362696 http://oos.soest.hawaii.edu/thredds/dodsC/paci...
223 3432967.640784 http://ecowatch.ncddc.noaa.gov/thredds/dodsC/n...
234 3077189.375788 http://ecowatch.ncddc.noaa.gov/thredds/dodsC/n...

In [27]:
url=[]
code=[]
for item in bad_data:
    url.append(item[0])
    code.append(item[1])

In [28]:
d={}
d['url']=url
d['code']=code
bad = pd.DataFrame(d)

In [29]:
bad.head()


Out[29]:
code url
0 404 http://www.neracoos.org/thredds/dodsC/UMO/DSG/...
1 404 http://www.neracoos.org/thredds/dodsC/UMO/DSG/...
2 404 http://www.neracoos.org/thredds/dodsC/UMO/DSG/...
3 404 http://www.neracoos.org/thredds/dodsC/UMO/DSG/...
4 404 http://www.neracoos.org/thredds/dodsC/UMO/DSG/...

In [70]:
cd /usgs/data2/notebook


/usgs/data2/notebook

In [71]:
bad.to_csv('bad.csv')

In [72]:
good_sorted['size']/=1e6

In [73]:
good_sorted.to_csv('good.csv')

In [74]:
['neracoos' in url for url in list(bad['url'].values)].count(1)


Out[74]:
23

In [75]:
['axiom' in url for url in list(bad['url'].values)].count(1)


Out[75]:
82

In [76]:
['caricoos' in url for url in list(bad['url'].values)].count(1)


Out[76]:
150

In [77]:
['secoora' in url for url in list(bad['url'].values)].count(1)


Out[77]:
21

In [78]:
pwd


Out[78]:
u'/usgs/data2/notebook'

In [79]:
!git add *.csv

In [80]:
!git commit -m 'updating csv'


[master 032029d] updating csv
 8 files changed, 160299 insertions(+)
 create mode 100644 bad.csv
 create mode 100644 gdp_texas_county_ccsm_a1f1.csv
 create mode 100644 gdp_texas_county_ccsm_a1fi.csv
 create mode 100644 gdp_texas_county_ccsm_b1.csv
 create mode 100644 gdp_texas_county_prism.csv
 create mode 100644 good.csv
 create mode 100644 sresa2_gfdl-cm2-1_1_Prcp.csv
 create mode 100644 sresb1_gfdl-cm2-1_1_Prcp.csv

In [81]:
git push


  File "<ipython-input-81-7e0a27daac38>", line 1
    git push
           ^
SyntaxError: invalid syntax

In [82]:
!git push


ssh: /home/usgs/anaconda/lib/libcrypto.so.1.0.0: no version information available (required by ssh)
Counting objects: 10, done.
Delta compression using up to 16 threads.
Compressing objects: 100% (9/9), done.
Writing objects: 100% (9/9), 680.67 KiB, done.
Total 9 (delta 3), reused 0 (delta 0)
To git@github.com:rsignell-usgs/notebook.git
   8598f3b..032029d  master -> master

In [ ]: