[GIS] Using WCS and Python to download from USDA Vegscape database

downloadowslibpythonwcswms

I wish to download Weekly NDVI rasters from the USDA Vegscape portal. It is possible to access the data using WMS and WCS methods. The USDA dev guide only contains example URLs on how to access the data, see below:

https://nassgeodata.gmu.edu/VegScapeService/weekly_ndvi_2012.cgi?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&coverage=Weekly_NDVI_30_2012.07.24_2012.07.30&CRS=epsg:102004&bbox=130783,2203171,153923,2217961&resx=250&resy=250&format=gtiff

I found the below section of code on gis.se but am unsure how to adapt it for my own purpose:

import os
from owslib.wcs import WebCoverageService as w  #owslib ->https://github.com/geopython/OWSLib
                                                #easy_install owslib

folder = 'c:\\path\\to\\folder\\'
wcs = w(url="http://dmcsee.org/cgi-bin/mapserv?map=/var/www/tmp/dmcsee_wms/dmcsee_wcs.map",version="1.0.0")
spi6 = w['SPI6']  #The layer I am interested in 

for timep in spi6.timepositions:
   output = wcs.getCoverage('SPI6',time=[timep],bbox=(10.0,31.995,48.005,50),format='GTiff',CRS='EPSG:4326', WIDTH=380, HEIGHT=120)
   f = open(os.path.join(folder,'spi6_'+str(timep)+'.tiff'),'wb')
   f.write(output.read())
   f.close()

Best Answer

You didn't say what you tried, or where you got stuck, but hopefully this answer will provide you enough hints to get further.

First off the original answer you referenced and included in your question had a syntax error; so as mentioned in comments:

for the ContentMetadata object you need w['SPI6'] not w.['SPI6']

Now looking at the service you want to access (https://nassgeodata.gmu.edu/VegScapeService/weekly_ndvi_2012.cgi?) when I try to call this using OWSLib I get the following error:

DEBUG:root:("bad handshake: Error([('SSL routines', 'ssl3_get_server_certificate', 'certificate verify failed')],)",)

Which means that there is a problem with the certificate used by the site. However, before we panic, if we look at the GetCapabilities response we can see that this service is proxying another:

<OnlineResource 
  xlink:type="simple" 
  xlink:href="http://ows9.csiss.gmu.edu/cgi-bin/daily_ndvi_2012?"/>

And we can use this URL in our Python like:

import owslib.wcs as wcs
nass_service = "http://ows9.csiss.gmu.edu/cgi-bin/daily_ndvi_2012?"
w = wcs.WebCoverageService(nass_service, version='1.0.0')

# Get a list of all the coverages in the service
print(list(w.contents))

Gives us:

'Daily_NDVI_2012.01.01', 'Daily_NDVI_2012.01.02', 'Daily_NDVI_2012.01.03', 'Daily_NDVI_2012.01.04', 'Daily_NDVI_2012.01.05', 'Daily_NDVI_2012.01.06', 'Daily_NDVI_2012.01.07', 'Daily_NDVI_2012.01.08'...

# Using one coverage (Daily_NDVI_2012.01.01) as an example
c = w['Daily_NDVI_2012.01.01']

# Take a look at the attributes of the coverage we could query
print(dir(c))

# Such as:
print("... c.c.boundingboxes ...")
for b in c.boundingboxes:
    print(b)

Gives:

... c.c.boundingboxes ...

{'nativeSrs': 'EPSG:4326', 'bbox': (-153.611110834706, 15.4492490834487, -33.1831608081618, 61.2417861047377)}

{'nativeSrs': 'epsg:102004', 'bbox': (-3987459.13480262, 168311.354467566, 4472862.72535489, 4177587.94694164)}

print("... c.timepositions ...")
for t in c.timepositions:
    print(t)

Gives:

... c.timepositions ...

print("... c.supportedFormats ...")
for f in c.supportedFormats:
  print(f)

Gives:

... c.supportedFormats ...

GTiff

So we can see there are no timepositions and therefore following an example where we loop over timepositions isn't going to be very useful.

except:

I realised that the https://service is actually proxying the daily service, but there is also also an accompanying weekly service we can call:

http://ows9.csiss.gmu.edu/cgi-bin/weekly_ndvi_2012?

However, the GetCapabilities response for this weekly service is also pointing at the daily service (like the https://service), so we have a problem because, as far as I know, it is not possible to get OWSLib to use the URL provided for wcs.WebCoverageService in the subsequent GetCoverage request, and instead the URLs from the GetCapabilites/DescribeCoverage responses are used...

I should add that generally this is the desired behaviour, but in this situation it is not helpful.

finally:

So if we want to get coverages from the weekly service, then we must write our own GetCoverage requests and in the end you'll need something like:

import tempfile
import os
import owslib.util
import owslib.wcs as wcs
import requests
import shutil

nass_daily_service = "http://ows9.csiss.gmu.edu/cgi-bin/daily_ndvi_2012?"
nass_weekly_service = "http://ows9.csiss.gmu.edu/cgi-bin/weekly_ndvi_2012?"

wcs_url = nass_weekly_service
w = wcs.WebCoverageService(wcs_url, version='1.0.0')
c = w['Weekly_NDVI_02_2012.01.10_2012.01.16']

resx = 1000
resy = 1000
rfmt = "GTiff"
rCRS = c.boundingboxes[1]['nativeSrs']
rBBOX = c.boundingboxes[1]['bbox']

tmpdir = tempfile.gettempdir()
fn = os.path.join(tmpdir,'test.tif')

try:
    gc = w.getCoverage(identifier=c.id, bbox=rBBOX, format=rfmt, crs=rCRS, width=resx, height=resy)
    f = open(fn, 'wb')
    f.write(gc.read())
    # This is what we would like to work
except owslib.util.ServiceException as se:
    #print(se)
    kvp_sep ="&"
    gc = str(nass_weekly_service) + "request=GetCoverage&service=WCS&version=1.0.0&Coverage=" + str(c.id) + "&bbox=" + str(rBBOX) + "&format=" +  str(rfmt) + "&width=" + str(resx) + "&height=" + str(resy) + "&crs=" + str(rCRS) + str(kvp_sep)

    r = requests.get(gc, stream=True)
    if r.status_code == 200:
        with open(fn, 'wb') as f:
            r.raw.decode_content = True
            shutil.copyfileobj(r.raw, f)
finally:
    f.close()
Related Question