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"
# 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
# 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
# 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()
# 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()
# 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
# 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()