NRCan and USGS Elevation products WCS extraction

In [1]:
import os
import gdal
import matplotlib.pyplot as plt
from matplotlib import image
import scipy.ndimage
import numpy as np
from numpy import*
import os
import folium
import rasterio
from rasterio.plot import plotting_extent
from ipyleaflet import Map, Marker, basemaps, basemap_to_tiles

USA_WCS = "https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer"
CANADA_WCS = "https://datacube.services.geo.ca/ows/elevation"
In [2]:
# Load OGC Services from NRCan and USGS

m = folium.Map(location=[65.0, -141.1],zoom_start=5)

canada = folium.raster_layers.WmsTileLayer(url = 'https://datacube.services.geo.ca/ows/elevation?', 
                                        layers = 'hrdsm-hillshade',
                                        fmt='image/png',
                                        version='1.3.0'
                                        )

usa = folium.raster_layers.WmsTileLayer(url = 'https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?',
                                       layers = '3DEPElevation:Hillshade Gray',
                                       fmt='image/png',
                                       version='1.3.0')
canada.layer_name = 'Canada'
usa.layer_name = 'USA'

usa.add_to(m)
canada.add_to(m)

folium.LayerControl().add_to(m)

m
Out[2]:
In [3]:
# Draw a polygon for extraction

xmin = -141
ymin = 64.6
xmax = -140.8
ymax = 64.8
aoi = """
{"type": "Polygon","coordinates": [[[%s,%s],[%s,%s],[%s,%s],[%s,%s],[%s,%s]]]}
""" % (xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin,xmin,ymin)

lat,lon = ((ymax+ymin)/2), ((xmax+xmin)/2)

m = folium.Map(location=[lat, lon],zoom_start=10)

usa.add_to(m)
canada.add_to(m)

folium.GeoJson(
    aoi,
    name='geojson'
).add_to(m)


folium.LayerControl().add_to(m)

m
Out[3]:
In [4]:
# Extract Canadian area from polygon (curl on WCS)

request = """%s?
SERVICE=WCS&
VERSION=1.1.1&
REQUEST=GetCoverage&
FORMAT=image/geotiff&
IDENTIFIER=dsm&
BOUNDINGBOX=%s,%s,%s,%s,urn:ogc:def:crs:EPSG::4326&
GRIDBASECRS=urn:ogc:def:crs:EPSG::4326&
GRIDCS=urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS&
GRIDTYPE=urn:ogc:def:method:WCS:1.1:2dSimpleGrid&
GRIDORIGIN=%s,%s&
GRIDOFFSETS=-0.00032835370435852,0.00032877413163684"""%(CANADA_WCS,ymin,xmin,ymax,xmax,ymax,xmin) 

request = "%s" % (request.replace('\n',''))
call = """curl -k -s "%s" --output Canada_tmp.tif"""%(request)

os.system(call)

call = "gdal_translate -ot Float32 Canada_tmp.tif Canada.tif"

os.system(call)

gtif = gdal.Open("Canada.tif")
georaster = gtif.ReadAsArray()

shape = georaster.shape

plt.imshow(georaster)
plt.title("Canada")
plt.show()
In [5]:
# Extract USA area from polygon (curl on WCS)

request = """%s?
SERVICE=WCS&
VERSION=1.0.0&
REQUEST=GetCoverage&
FORMAT=GeoTIFF&
COVERAGE=DEP3Elevation&
BBOX=%s,%s,%s,%s&
CRS=EPSG:4326&
RESPONSE_CRS=EPSG:4326&
WIDTH=%s&
HEIGHT=%s"""%(USA_WCS,xmin,ymin,xmax,ymax,shape[1],shape[0]) 

request = "%s" % (request.replace('\n',''))
call = """curl -k -s "%s" --output USA.tif"""%(request)

os.system(call)

gtif = gdal.Open("USA.tif")
georaster = gtif.ReadAsArray()

plt.imshow(georaster)
plt.title("USA")
plt.show()
In [6]:
# Plot differences between two raster

with rasterio.open("Canada.tif") as canada_dem:
    canada_dem_im = canada_dem.read(1,masked=True)
    bounds = plotting_extent(canada_dem)

with rasterio.open("USA.tif") as usa_dem:
    usa_dem_im = usa_dem.read(1,masked=True)

# Are the bounds the same?
print("Is the spatial extent the same?", 
      canada_dem.bounds == usa_dem.bounds)

# Is the resolution the same ??
print("Is the resolution the same?", 
      canada_dem.res == usa_dem.res)

diff = canada_dem_im - usa_dem_im

plt.imshow(diff)
plt.title("Difference Canada/USA")
plt.colorbar()
plt.show()

plt.hist(diff)
plt.show()
D:\apps\conda\envs\datacube\lib\site-packages\ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in greater
  after removing the cwd from sys.path.
D:\apps\conda\envs\datacube\lib\site-packages\ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in less
  after removing the cwd from sys.path.
Is the spatial extent the same? False
Is the resolution the same? False
In [7]:
# Merge DTMs, draw a line, extract differences

# Generate some data
call = "gdalbuildvrt Canada-USA-merged.vrt USA.tif Canada.tif"

os.system(call)

gtif = gdal.Open("Canada-USA-merged.vrt")
georaster = gtif.ReadAsArray()
i = isnan(georaster)
georaster[i] = 0

canada = gdal.Open("Canada.tif")
canada_grid = canada.ReadAsArray()
i = isnan(canada_grid)
canada_grid[i] = 0

usa = gdal.Open("USA.tif")
usa_grid = usa.ReadAsArray()


# Extract the line
x0, y0 = 225, 150
x1, y1 = 55,400
num = 101
x, y = np.linspace(x0,x1,num), np.linspace(y0,y1,num)

# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(georaster, np.vstack((x,y)))
canada = scipy.ndimage.map_coordinates(canada_grid, np.vstack((x,y)))
usa = scipy.ndimage.map_coordinates(usa_grid, np.vstack((x,y)))

# Plot...

fig, axes= plt.subplots(nrows=3, figsize=(20,30))

axes[0].imshow(georaster, interpolation='none')
axes[0].plot([x0, x1], [y0, y1], 'ro-')
axes[0].axis('image')
axes[0].set_title('Canada/USA merged')

l1 = axes[1].plot(canada, 'r',linewidth=3, label='Canada')
l2 = axes[1].plot(usa,'b', linewidth=3, label='USA')
l3 = axes[1].plot(zi, '--g', linewidth=3, label='Canada-USA merged')
axes[1].legend()

l1 = axes[2].plot(canada, 'r',linewidth=3, label='Canada')
l2 = axes[2].plot(usa,'b', linewidth=3, label='USA')
l3 =axes[2].plot(zi, '--g', linewidth=5, label='Canada-USA merged')
axes[2].set_xlim([30,60])
axes[2].set_ylim([400,1000])
axes[2].legend()

plt.show()
In [ ]: