from shapely.geometry import LineString
# Dilating a line
line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
dilated = line.buffer(0.5)
eroded = dilated.buffer(-0.3)
# Demonstate Python Geo Interface
print(line.__geo_interface__)
The data was originally sourced from the HURDAT2 dataset from AOML/NOAA, and the Python lists are from the cartopy documentation.
from shapely.geometry import LineString
lons = [-75.1, -75.7, -76.2, -76.5, -76.9, -77.7, -78.4, -79.0,
-79.6, -80.1, -80.3, -81.3, -82.0, -82.6, -83.3, -84.0,
-84.7, -85.3, -85.9, -86.7, -87.7, -88.6, -89.2, -89.6,
-89.6, -89.6, -89.6, -89.6, -89.1, -88.6, -88.0, -87.0,
-85.3, -82.9]
lats = [23.1, 23.4, 23.8, 24.5, 25.4, 26.0, 26.1, 26.2, 26.2, 26.0,
25.9, 25.4, 25.1, 24.9, 24.6, 24.4, 24.4, 24.5, 24.8, 25.2,
25.7, 26.3, 27.2, 28.2, 29.3, 29.5, 30.2, 31.1, 32.6, 34.1,
35.6, 37.0, 38.6, 40.1]
# Turn the lons and lats into a shapely LineString
katrina_track = LineString(zip(lons, lats))
Buffer the linestring by two degrees (doesn't really make sense!). This should be about 200kms, but as we'll see, it's not quite accurate... Why not?
katrina_buffer = katrina_track.buffer(2)
What if we reproject the lon/lats to a projection that preserves distances better?
We could use EPSG:32616
(UTM Zone 16), which covers where Katrina meets New Orleans, but we're probably better off using a custom proj4
string based on a Lambert Conformal Conic projection. Why?
from pyproj import Proj, transform
from fiona.crs import from_string
# Create custom proj4 string
proj = from_string('+ellps=WGS84 +proj=lcc +lon_0=-96.0 +lat_0=39.0 '
'+x_0=0.0 +y_0=0.0 +lat_1=33 +lat_2=45 +no_defs')
# Create source and destination Proj objects (source is WGS84 lons/lats)
src = Proj(init='epsg:4326')
dst = Proj(proj)
# Create a LineString from the transformed points
proj_track = LineString(zip(*transform(src, dst, lons, lats)))
# Buffer the LineString by 200 km (multiply by 1000 to work in meters)
proj_buffer = proj_track.buffer(200*1000)
zip
is your friend! Use it with tuple unpacking to change between sequences of (x, y)
pairs and seperate x
and y
sequences.
pts = [(0, 0), (1, 0), (1, 1), (2, 1), (2, 2)]
x, y = zip(*pts)
print x, y
Also, instead of calling f(x, y)
, you can just use
# Skip this slide, not really needed for demo, just need Python function
# Simple function to add 0.5 to each coordinate
def f(x, y):
new_x = [i + 0.5 for i in x]
new_y = [j + 0.5 for j in y]
return new_x, new_y
f(*zip(*pts)) # Function f adds 0.5 to each coordinate
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.LambertConformal())
ax.stock_img() # Add background image (slow)
ax.coastlines(resolution='110m')
ax.add_geometries([katrina_buffer], ccrs.PlateCarree(), facecolor='blue', alpha=0.5)
ax.add_geometries([katrina_track], ccrs.PlateCarree(), facecolor='none')
# Let's add the projected buffer for comparison
ax.add_geometries([proj_buffer], ccrs.LambertConformal(), facecolor='green', alpha=0.5)
ax.set_extent([-125, -66.5, 20, 50], ccrs.PlateCarree())
ax.gridlines()
plt.show()
Which shapely
geometry method could we use to find where the tracks differ?
What makes geospatial raster datasets different from other raster files is that their pixels map to regions of the Earth. In this case, we have a raster image which maps to Midtown Manhattan.
import matplotlib.image as mpimg
import os
# Read in a regular PNG image of Manhattan
png_file = os.path.join('..', 'data', 'manhattan.png')
img = mpimg.imread(png_file)
# Take a quick look at the shape
print(img.shape)
# Specify the affine transformation
from affine import Affine # You might not have this library installed
A = Affine(0.999948245999997, 0.0, 583057.357, 0.0, -0.999948245999997, 4516255.36)
# Compute the upper left and lower right corners
ul = (0, 0)
lr = img.shape[:2][::-1]
print("Upper left: " + str(A * ul))
print("Lower right: " + str(A * lr))
Here, the coordinate reference system is EPSG:26918
(NAD83 / UTM Zone 18N), and the affine transformation matrix is given (in later examples, we'll get this information directly from the input raster datasets).
# Upper left and bottom right corners in UTM coords
left, top = A * ul
right, bottom = A * lr
# Plot showing original PNG image (axes correspond to rows and cols) on left
# and 'transformed' PNG (axes correspond to UTM coords) on the right
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))
ax1.imshow(img)
ax1.set_title("PNG with row, col bounds")
ax2.imshow(img, extent=(left, right, bottom, top), aspect="equal")
ax2.set_title("PNG with correct bounds")
plt.show()
[Plotting Great Circles in Python
](../exercises/Plotting Great Circles in Python.ipynb)