Visualizing DEM differencing data

Objective: create an interactive map that has a dem or hillshade as a basemap and allows the user to use a scrollbar to view DEMs and DEM difference geotifs

The following python packages are required:

Note: this is a work in progress. The creation of an interactive map with a DEM overlay part has been completed, but other sections are still in development. Some code blocks are marked "do not run" because they are not needed for the DEM overlay

import python libraries


In [1]:
#do not run
import sys
#update path until georaster is installed with the make file
sys.path.insert(0,'/Users/jessica/Classes/Geohackweek2016/iceflow/georaster')

In [2]:
#do not run
import geoutils
import gdal
import pandas
from matplotlib import pyplot as plt


//anaconda/envs/py27/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

In [3]:
import mpl_toolkits.basemap
from ipyleaflet import (Map,
    Marker,
    TileLayer, ImageOverlay,
    Polyline, Polygon, Rectangle, Circle, CircleMarker,
    GeoJSON,
    DrawControl
)

Set path to input raster and read in for visualizing


In [5]:
#import georaster as raster  ##sometimes needs to be from georaster import georaster...
#rast_img = raster.SingleBandRaster('../iceflow_project/20141022_0518_1020010036518E00_102001003525D400-DEM_32m.tif')

In [8]:
#do not run
%matplotlib inline

Set up working environment/variables


In [25]:
#do not run - a work in progress to not have the corners hard coded in
ds = gdal.Open('../iceflow/data/WV_Ngozumpa/coreg32_mos-tile-0.tif')
ulx, xres, xskew, uly, yskew, yres = ds.GetGeoTransform()

llx = ulx
lly = uly - (ds.RasterYSize * yres)
urx = ulx + (ds.RasterXSize * xres)
ury = uly

In [46]:
#do not run - a work in progress to not have the corners hard coded in
import subprocess
cmd ="gdalinfo ./data/WV_Ngozumpa/coreg32_mos-tile-0.tif"
print subprocess.check_output(cmd,shell=True)
print os.system('pwd')


Driver: GTiff/GeoTIFF
Files: ./data/WV_Ngozumpa/coreg32_mos-tile-0.tif
Size is 608, 1958
Coordinate System is:
PROJCS["WGS 84 / UTM zone 45N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32645"]]
Origin = (460048.000000000000000,3146768.000000000000000)
Pixel Size = (32.000000000000000,-32.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  460048.000, 3146768.000) ( 86d35'31.09"E, 28d26'48.54"N)
Lower Left  (  460048.000, 3084112.000) ( 86d35'38.80"E, 27d52'52.43"N)
Upper Right (  479504.000, 3146768.000) ( 86d47'26.42"E, 28d26'50.17"N)
Lower Right (  479504.000, 3084112.000) ( 86d47'30.37"E, 27d52'54.02"N)
Center      (  469776.000, 3115440.000) ( 86d41'31.69"E, 28d 9'51.44"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=0

0

In [55]:
center = [28.1, 86.6] #center coordinates, lat then lon, for the base interactive map
zoom = 10
DEM_url = 'http://i.imgur.com/5DK99mp.jpg'
bound_SW = [27.881231, 86.594111]
bound_NE = [28.447269, 86.790672]

Create map, then add an overlay


In [56]:
map = Map(center=center, zoom=zoom)
#map

Important notes about ImageOverlay: the url MUST point to something on the web. It does not work with local addresses (relative or absolute). Additionally, it does not work with services that require a login (e.g. google drive -> with a shareable link to the drive combined with a path or a shareable link to the image directly, github). It was discovered that putting an image on imgur and then using that link works!

Attempts at an actual dem dataset/image:
- https://drive.google.com/drive/folders/0B5c3UTO8DDZwamlQMzNRRWd3TEU/HMA/incoming/COREG/coreg_32m/20151018_0514_1050010001B64900_1050010001B64A00-DEM_32m_fig.png'

Example of a working line of code to add an overlay: DEM_overlay = ImageOverlay(url='https://s3.amazonaws.com/images.seroundtable.com/t-google-star-1300195846.jpg', bounds=[[27.8, 86.9],[28, 87.3]])


In [58]:
map.remove_layer(DEM_overlay)
DEM_overlay = ImageOverlay(url=DEM_url, bounds=[bound_SW, bound_NE],\
                          opacity=0.7)
map.add_layer(DEM_overlay)
map

In [ ]: