GeoHackWeek 2019 -- vector tutorial
Friedrich Knuth, University of Washington. 2019-9-9
import cartopy.crs as ccrs
import holoviews as hv
import geoviews as gv
from geoviews import opts
from holoviews.streams import PointDraw
import panel as pn
import warnings
warnings.filterwarnings('ignore')
hv.extension('bokeh')
# On pangeo, uncomment this line to install this package not present in the pangeo conda environment
# !conda install -y -c conda-forge contextily
OpenTopoMap = 'https://tile.opentopomap.org/{Z}/{X}/{Y}.png'
OpenStreetMap = 'http://tile.openstreetmap.org/{Z}/{X}/{Y}.png'
GoogleHybrid = 'https://mt1.google.com/vt/lyrs=y&x={X}&y={Y}&z={Z}'
GoogleSatellite = 'https://mt1.google.com/vt/lyrs=s&x={X}&y={Y}&z={Z}'
GoogleRoad = 'https://mt1.google.com/vt/lyrs=m&x={X}&y={Y}&z={Z}'
GoogleTerrain = 'http://mt0.google.com/vt/lyrs=p&hl=en&x={X}&y={Y}&z={Z}'
GoogleTerrainOnly = 'http://mt0.google.com/vt/lyrs=t&hl=en&x={X}&y={Y}&z={Z}'
ESRIImagery = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{Z}/{Y}/{X}.jpg'
Wikimedia = 'https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png'
url = GoogleSatellite
lon = -122.3118
lat = 47.6534
delta = 0.001
extents = (lon-delta, lat-delta, lon+delta, lat+delta)
tiles = gv.WMTS(url, extents=extents)
location = gv.Points([], vdims="vertices_label")
point_stream = PointDraw(source=location)
base_map = (tiles * location).opts(opts.Points(width=500,
height=500,
size=5,
color='blue',
tools=["hover"]))
app = pn.panel(base_map)
app
df = gv.operation.project_points(point_stream.element).dframe()
df
x | y | vertices_label | |
---|---|---|---|
0 | -1.361570e+07 | 6.049399e+06 | None |
1 | -1.361570e+07 | 6.049369e+06 | None |
2 | -1.361567e+07 | 6.049369e+06 | None |
3 | -1.361567e+07 | 6.049397e+06 | None |
import contextily
from contextily.tile import bounds2raster
url = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}"
output_file_name = 'output.tif'
w, s, e, n = (lon-delta, lat-delta, lon+delta, lat+delta)
img = bounds2raster(
w,
s,
e,
n,
output_file_name,
zoom="auto",
url=url,
ll=True,
wait=0,
max_retries=100)
! gdalinfo output.tif
Driver: GTiff/GeoTIFF Files: output.tif Size is 768, 1536 Coordinate System is: PROJCRS["WGS 84 / Pseudo-Mercator", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["Popular Visualisation Pseudo-Mercator", METHOD["Popular Visualisation Pseudo Mercator", ID["EPSG",1024]], PARAMETER["Latitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",0, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["False easting",0, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",0, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["easting (X)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (Y)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["unknown"], AREA["World - 85°S to 85°N"], BBOX[-85.06,-180,85.06,180]], ID["EPSG",3857]] Data axis to CRS axis mapping: 1,2 Origin = (-13615804.434757798910141,6049608.752921344712377) Pixel Size = (0.298582141738734,-0.298582141739341) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left (-13615804.435, 6049608.753) (122d18'46.27"W, 47d39'17.10"N) Lower Left (-13615804.435, 6049150.131) (122d18'46.27"W, 47d39' 7.11"N) Upper Right (-13615575.124, 6049608.753) (122d18'38.85"W, 47d39'17.10"N) Lower Right (-13615575.124, 6049150.131) (122d18'38.85"W, 47d39' 7.11"N) Center (-13615689.779, 6049379.442) (122d18'42.56"W, 47d39'12.11"N) Band 1 Block=768x3 Type=Byte, ColorInterp=Red Band 2 Block=768x3 Type=Byte, ColorInterp=Green Band 3 Block=768x3 Type=Byte, ColorInterp=Blue
import rasterio
import matplotlib.pyplot as plt
from rasterio.plot import show
source = rasterio.open(output_file_name)
fig, ax = plt.subplots(figsize=(10, 10))
show(source.read(),
transform=source.transform,
ax=ax)
ax.ticklabel_format(useOffset=False, style='plain', axis='both')
from shapely.geometry import Point, Polygon
import geopandas as gpd
# this is the projection the tiles are returned in
epsg_code = '3857'
geometry = [Point(xy) for xy in zip(df['x'], df['y'])]
gdf = gpd.GeoDataFrame(df,geometry=geometry, crs={'init':'epsg:'+epsg_code})
fig, ax = plt.subplots(figsize=(10, 10))
show(source.read(),
transform=source.transform,
ax=ax)
ax.ticklabel_format(useOffset=False, style='plain', axis='both')
gdf.plot(ax=ax, color='blue')
<matplotlib.axes._subplots.AxesSubplot at 0x7fd4d50f1050>
vertices = []
for i in range(len(df)):
vertex = (df['x'][i], df['y'][i])
vertices.append(vertex)
polygon = Polygon(vertices)
polygon_gdf = gpd.GeoDataFrame(gpd.GeoSeries(polygon),
columns=['geometry'],
crs={'init':'epsg:'+epsg_code})
fig, ax = plt.subplots(figsize=(10, 10))
show(source.read(),
transform=source.transform,
ax=ax)
ax.ticklabel_format(useOffset=False, style='plain', axis='both')
gdf.plot(ax=ax, color='blue')
polygon_gdf.plot(ax=ax)
<matplotlib.axes._subplots.AxesSubplot at 0x7fd4d40ad550>
# ! pip install utm
import utm
url = 'https://mt1.google.com/vt/lyrs=s&x={X}&y={Y}&z={Z}'
u = utm.from_latlon(lat,lon)
utm_lon = u[0]
utm_lat = u[1]
utm_zone = u[2]
utm_zone_code = u[3]
dx = dy = 100
extents = utm_lon-dx, utm_lat-dy, utm_lon+dx, utm_lat+dy
tiles = gv.WMTS(url, extents=extents, crs=ccrs.UTM(utm_zone))
location = gv.Points([], vdims="location", crs=ccrs.UTM(utm_zone))
point_stream = PointDraw(source=location)
base_map = (tiles * location).opts(gv.opts.Points(width=500,
height=500,
size=5,
color='blue',
tools=["hover"]))
app = pn.panel(base_map)
app
projected = gv.operation.project_points(point_stream.element, projection=ccrs.UTM(utm_zone))
projected.dframe()
x | y | location | |
---|---|---|---|
0 | 551668.638822 | 5.278018e+06 | None |
1 | 551669.328656 | 5.277998e+06 | None |
2 | 551690.722934 | 5.277998e+06 | None |
3 | 551690.574091 | 5.278015e+06 | None |
import math
# function to get utm epsg code from wgs84 lon lat values
def wgs_lon_lat_to_utm_epsg_code(lon, lat):
utm_band = str((math.floor((lon + 180) / 6 ) % 60) + 1)
if len(utm_band) == 1:
utm_band = '0'+utm_band
if lat >= 0:
epsg_code = '326' + utm_band
else:
epsg_code = '327' + utm_band
return epsg_code
epsg_code = wgs_lon_lat_to_utm_epsg_code(lon, lat)
reprojected_output_file_name = 'output_EPSG_'+epsg_code+'.tif'
! rm {reprojected_output_file_name}
# change the projection using gdal
! gdalwarp -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER -r cubic -t_srs EPSG:{epsg_code} {output_file_name} {reprojected_output_file_name}
Creating output file that is 783P x 1542L. Processing output.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.
src = rasterio.open(reprojected_output_file_name)
fig, ax = plt.subplots(figsize=(10, 10))
show(src.read(),
transform=src.transform,
ax=ax)
ax.ticklabel_format(useOffset=False, style='plain', axis='both')
Can you reproject and plot the vertices and ploygon on top of the UTM projected raster?
Hint: Geopandas has a built in function to do this. You will need to provide the correct epsg code, which you can get with the wgs_lon_lat_to_utm_epsg_code(lon, lat) function defined in this notebook.
See further examples at https://holoviz.org/tutorial/index.html