Week 4B: Working with Raster Data¶
Sep 29, 2021
Housekeeping¶
- Homework #2 due today
- Homework #3 posted today, due in two weeks (10/13)
- Part 1: Working with eviction data in Philadelphia
- Part 2: Exploring the NDVI in Philadelphia
- Come to office hours!
- Tuesday evening (Nick) 6-8pm
- Monday afternoon (Stella) 12:30 pm - 2 pm
- Sign up for time slots on Canvas calendar
Follow along in lecture on Binder!¶
All lecture slides will work on Binder, data is available to be loaded, etc. Binder has an exact copy of the file structure in the GitHub repo.
https://github.com/MUSA-550-Fall-2021/week-4
Last time: A set of coordinated visualization libraries in Python¶
hvplot¶
- Quickly generate interactive plots from your data
- Seamlessly handles pandas and geopandas data
- Relies on Holoviews and Geoviews under the hood
- Fancy, magic widgets for grouping charts by a third column
- Great built-in interaction toolbar powered by bokeh: box select, lasso, pan, zoom, etc...
An interface just like the pandas plot() function, but much more useful.
In [8]:
# Our usual imports
import pandas as pd
import geopandas as gpd
import numpy as np
import altair as alt
from matplotlib import pyplot as plt
In [9]:
# This will add the .hvplot() function to your DataFrame!
import hvplot.pandas
# Import holoviews too
import holoviews as hv
Last time: hvplot has great support for geopandas¶
Example: Interactive choropleth of median property assessments in Philadelphia
In [10]:
# Load the data
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2021/week-3/main/data/opa_residential.csv"
data = pd.read_csv(url)
# Create the Point() objects
data['geometry'] = gpd.points_from_xy(data['lng'], data['lat'])
# Create the GeoDataFrame
data = gpd.GeoDataFrame(data, geometry='geometry', crs="EPSG:4326")
In [11]:
# load the Zillow data from GitHub
url = "https://raw.githubusercontent.com/MUSA-550-Fall-2021/week-3/main/data/zillow_neighborhoods.geojson"
zillow = gpd.read_file(url)
In [12]:
# Important: Make sure the CRS match
data = data.to_crs(zillow.crs)
# perform the spatial join
data = gpd.sjoin(data, zillow, op='within', how='left')
In [13]:
# Calculate the median market value per Zillow neighborhood
median_values = data.groupby("ZillowName", as_index=False)["market_value"].median()
# Merge median values with the Zillow geometries
median_values = zillow.merge(median_values, on="ZillowName")
print(type(median_values))
In [14]:
median_values.head()
Out[14]:
Take advantage of GeoViews¶
Let's add a tile source underneath the choropleth map
In [15]:
import geoviews as gv
import geoviews.tile_sources as gvts
In [16]:
%%opts WMTS [width=800, height=800, xaxis=None, yaxis=None]
choro = median_values.hvplot(c='market_value',
width=500,
height=400,
alpha=0.5,
geo=True,
cmap='viridis',
hover_cols=['ZillowName'])
gvts.ESRI * choro
Out[16]:
hvplot also works with raster data¶
- So far we've been working mainly with vector data using geopandas: lines, points, polygons
- The basemaps we've been adding are an exampe of raster data
- Raster data:
- Gridded or pixelated
- Maps easily to an array — it's the representation for digital images
Continuous raster examples¶
- Multispectral satellite imagery, e.g., the Landsat satellite
- Digital Elevation Models (DEMs)
- Maps of canopy height derived from LiDAR data.
Categorical raster examples¶
- Land cover or land-use maps.
- Tree height maps classified as short, medium, and tall trees.
- Snowcover masks (binary snow or no snow)
Important attributes of raster data¶
1. The coordinate reference system¶
4. Affine transformation¶
Pixel space --> real space
- A mapping from pixel space (rows and columns of the image) to the x and y coordinates defined by the CRS
- Typically a six parameter matrix defining the origin, pixel resolution, and rotation of the raster image
Multi-band raster data¶
Each band measures the light reflected from different parts of the electromagnetic spectrum.
Color images are multi-band rasters!¶
The raster format: GeoTIFF¶
- A standard
.tifimage format with additional spatial information embedded in the file, which includes metadata for:- Geotransform, e.g., extent, resolution
- Coordinate Reference System (CRS)
- Values that represent missing data (
NoDataValue)
Tools for raster data¶
- Lots of different options: really just working with arrays
- One of the first options: Geospatial Data Abstraction Library (GDAL)
- Low level and not really user-friendly
- A more recent, much more Pythonic package:
rasterio
We'll use rasterio for the majority of our raster analysis
Getting started with rasterio¶
Analysis is built around the "open()" command which returns a "dataset" or "file handle"
In [17]:
import rasterio as rio
In [18]:
# Open the file and get a file "handle"
landsat = rio.open("./data/landsat8_philly.tif")
landsat
Out[18]:
Let's check out the metadata¶
In [19]:
# The CRS
landsat.crs
Out[19]:
In [20]:
# The bounds
landsat.bounds
Out[20]:
In [21]:
# The number of bands available
landsat.count
Out[21]:
In [22]:
# The band numbers that are available
landsat.indexes
Out[22]:
In [23]:
# Number of pixels in the x and y directions
landsat.shape
Out[23]:
In [24]:
# The 6 parameters that map from pixel to real space
landsat.transform
Out[24]:
In [25]:
# All of the meta data
landsat.meta
Out[25]:
Reading data from a file¶
Tell the file which band to read, starting with 1
In [26]:
# This is just a numpy array
data = landsat.read(1)
data
Out[26]:
We can plot it with matplotlib's imshow
In [27]:
fig, ax = plt.subplots(figsize=(10,10))
img = ax.imshow(data)
plt.colorbar(img)
Out[27]:
Two improvements¶
- A log-scale colorbar
- Add the axes extent
Adding the axes extent¶
Note that imshow needs a bounding box of the form: (left, right, bottom, top).
See the documentation for imshow if you forget the right ordering...(I often do!)
Using a log-scale color bar¶
Matplotlib has a built in log normalization...
In [28]:
import matplotlib.colors as mcolors
In [29]:
# Initialize
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the image
img = ax.imshow(
data,
norm=mcolors.LogNorm(), # Use a log colorbar scale
extent=[ # Set the extent of the images
landsat.bounds.left,
landsat.bounds.right,
landsat.bounds.bottom,
landsat.bounds.top,
],
)
# Add a colorbar
plt.colorbar(img)
Out[29]:
Overlaying vector geometries on raster plots¶
Two requirements:
- The CRS of the two datasets must match
- The extent of the imshow plot must be set properly
Let's add the city limits¶
In [30]:
city_limits = gpd.read_file("./data/City_Limits.geojson")
In [31]:
# Convert to the correct CRS!
print(landsat.crs.data)
In [32]:
city_limits = city_limits.to_crs(landsat.crs.data['init'])
In [33]:
# Initialize
fig, ax = plt.subplots(figsize=(10, 10))
# The extent of the data
landsat_extent = [
landsat.bounds.left,
landsat.bounds.right,
landsat.bounds.bottom,
landsat.bounds.top,
]
# Plot!
img = ax.imshow(data,
norm=mcolors.LogNorm(), #NEW
extent=landsat_extent) #NEW
# Add the city limits
city_limits.plot(ax=ax, facecolor="none", edgecolor="white")
# Add a colorbar and turn off axis lines
plt.colorbar(img)
ax.set_axis_off()
What band did we plot??¶
Landsat 8 has 11 bands spanning the electromagnetic spectrum from visible to infrared.
How about a typical RGB color image?¶
First let's read the red, blue, and green bands (bands 4, 3, 2):
Note¶
The .indexes attributes tells us the bands available to be read.
Important: they are not zero-indexed!
In [34]:
landsat.indexes
Out[34]:
In [35]:
rgb_data = landsat.read([4, 3, 2])
In [36]:
rgb_data.shape
Out[36]:
Note the data has the shape: (bands, height, width)
In [37]:
data[2]
Out[37]:
In [38]:
fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(12, 6))
cmaps = ["Reds", "Greens", "Blues"]
for i in [0, 1, 2]:
# This subplot axes
ax = axs[i]
# Plot
ax.imshow(rgb_data[i], norm=mcolors.LogNorm(), extent=landsat_extent, cmap=cmaps[i])
city_limits.plot(ax=ax, facecolor="none", edgecolor="black")
# Format
ax.set_axis_off()
ax.set_title(cmaps[i][:-1])
Side note: using subplots in matplotlib¶
You can specify the subplot layout via the nrows and ncols keywords to plt.subplots. If you have more than one row or column, the function returns:
- The figure
- The list of axes.
In [39]:
fig, axs = plt.subplots(nrows=1, ncols=3)
Note¶
When nrows or ncols > 1, I usually name the returned axes variable "axs" vs. the usual case of just "ax". It's useful for remembering we got more than one Axes object back!
In [40]:
# This has length 3 for each of the 3 columns!
axs
Out[40]:
In [41]:
# Left axes
axs[0]
Out[41]:
In [42]:
# Middle axes
axs[1]
Out[42]:
In [43]:
# Right axes
axs[2]
Out[43]:
Use earthpy to plot the combined color image¶
A useful utility package to perform the proper re-scaling of pixel values
In [44]:
import earthpy.plot as ep
In [45]:
# Initialize
fig, ax = plt.subplots(figsize=(10,10))
# Plot the RGB bands
ep.plot_rgb(rgb_data, rgb=(0, 1, 2), ax=ax);
We can "stretch" the data to improve the brightness¶
In [46]:
# Initialize
fig, ax = plt.subplots(figsize=(10,10))
# Plot the RGB bands
ep.plot_rgb(rgb_data, rgb=(0, 1, 2), ax=ax, stretch=True);
Conclusion: Working with multi-band data is a bit messy¶
Can we do better?
Yes! HoloViz to the rescue...¶
xarray¶
- A fancier version of numpy arrays
- Designed for gridded data with multiple dimensions
- For raster data, those dimensions are: bands, latitude, longitude, pixel values
Let's try it out...¶
In [47]:
import xarray as xr
In [48]:
xr.open_rasterio?
In [49]:
ds = xr.open_rasterio("./data/landsat8_philly.tif")
ds
Out[49]:
In [50]:
type(ds)
Out[50]:
Now the magic begins: more hvplot¶
Still experimental, but very promising — interactive plots with widgets for each band!
In [51]:
import hvplot.xarray
In [52]:
img = ds.hvplot.image(width=700, height=500, cmap='viridis')
img
Out[52]:
More magic widgets!¶
Two key raster concepts¶
- Using vector + raster data
- Cropping or masking raster data based on vector polygons
- Raster math and zonal statistics
- Combining multiple raster data sets
- Calculating statistics within certain vector polygons
Cropping: let's trim to just the city limits¶
Use rasterio's builtin mask function:
In [53]:
from rasterio.mask import mask
In [54]:
masked, mask_transform = mask(
dataset=landsat,
shapes=city_limits.geometry,
crop=True, # remove pixels not within boundary
all_touched=True, # get all pixels that touch the boudnary
filled=False, # do not fill cropped pixels with a default value
)
In [55]:
masked
Out[55]:
In [56]:
masked.shape
Out[56]:
In [57]:
# Initialize
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the first band
ax.imshow(masked[0], cmap="viridis", extent=landsat_extent)
# Format and add the city limits
city_limits.boundary.plot(ax=ax, color="gray", linewidth=4)
ax.set_axis_off()
Writing GeoTIFF files¶
Let's save our cropped raster data. To save raster data, we need both the array of the data and the proper meta-data.
In [58]:
out_meta = landsat.meta
out_meta.update(
{"height": masked.shape[1], "width": masked.shape[2], "transform": mask_transform}
)
print(out_meta)
# write small image to local Geotiff file
with rio.open("data/cropped_landsat.tif", "w", **out_meta) as dst:
dst.write(masked)
# Can't access
In [59]:
landsat
Out[59]:
Note: we used a context manager above (the "with" statement) to handle the opening and subsequent closing of our new file.
For more info, see this tutorial.
Now let's do some raster math¶
Often called "map algebra"...
The Normalized Difference Vegetation Index¶
- A normalized index of greenness ranging from -1 to 1, calculated from the red and NIR bands.
- Provides a measure of the amount of live green vegetation in an area
Formula:
NDVI = (NIR - Red) / (NIR + Red)
Healthy vegetation reflects NIR and absorbs visible, so the NDVI for green vegatation
Get the data for the red and NIR bands¶
NIR is band 5 and red is band 4
In [60]:
masked.shape
Out[60]:
In [61]:
# Note that the indexing here is zero-based, e.g., band 1 is index 0
red = masked[3]
nir = masked[4]
In [62]:
def calculate_NDVI(nir, red):
"""
Calculate the NDVI from the NIR and red landsat bands
"""
# Convert to floats
nir = nir.astype(float)
red = red.astype(float)
# Get valid entries
check = np.logical_and( red.mask == False, nir.mask == False )
# Where the check is True, return the NDVI, else return NaN
ndvi = np.where(check, (nir - red ) / ( nir + red ), np.nan )
return ndvi
In [63]:
NDVI = calculate_NDVI(nir, red)
In [64]:
fig, ax = plt.subplots(figsize=(10,10))
# Plot NDVI
img = ax.imshow(NDVI, extent=landsat_extent)
# Format and plot city limits
city_limits.plot(ax=ax, edgecolor='gray', facecolor='none', linewidth=4)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in Philadelphia", fontsize=18);
Let's overlay Philly's parks¶
In [65]:
# Read in the parks dataset
parks = gpd.read_file('./data/parks.geojson')
In [66]:
# Print out the CRS
parks.crs
Out[66]:
In [67]:
landsat.crs
Out[67]:
In [68]:
landsat.crs.data['init']
Out[68]:
In [69]:
# Conver to landsat CRS
parks = parks.to_crs(landsat.crs.data['init'])
#parks = parks.to_crs(epsg=32618)
In [70]:
fig, ax = plt.subplots(figsize=(10,10))
# Plot NDVI
img = ax.imshow(NDVI, extent=landsat_extent)
# NEW: add the parks
parks.plot(ax=ax, edgecolor='crimson', facecolor='none', linewidth=2)
# Format and plot city limits
city_limits.plot(ax=ax, edgecolor='gray', facecolor='none', linewidth=4)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI vs. Parks in Philadelphia", fontsize=18);
It looks like it worked pretty well!
How about calculating the median NDVI within the park geometries?¶
This is called zonal statistics. We can use the rasterstats package to do this...
In [71]:
from rasterstats import zonal_stats
In [72]:
stats = zonal_stats(parks, NDVI, affine=landsat.transform, stats=['mean', 'median'])
The zonal_stats function¶
Zonal statistics of raster values aggregated to vector geometries.
- First argument: vector data
- Second argument: raster data to aggregated
- Also need to pass the affine transform and the names of the stats to compute
In [73]:
stats
Out[73]:
In [74]:
len(stats)
Out[74]:
In [75]:
len(parks)
Out[75]:
In [76]:
# Store the median value in the parks data frame
parks['median_NDVI'] = [s['median'] for s in stats]
In [77]:
parks.head()
Out[77]:
Make a quick histogram of the median values¶
They are all positive, indicating an abundance of green vegetation...
In [78]:
# Initialize
fig, ax = plt.subplots(figsize=(8,6))
# Plot a quick histogram
ax.hist(parks['median_NDVI'], bins='auto')
ax.axvline(x=0, c='k', lw=2)
# Format
ax.set_xlabel("Median NDVI", fontsize=18)
ax.set_ylabel("Number of Parks", fontsize=18);
Let's make a choropleth, too¶
In [79]:
# Initialize
fig, ax = plt.subplots(figsize=(10,10))
# Plot the city limits
city_limits.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=4)
# Plot the median NDVI
parks.plot(column='median_NDVI', legend=True, ax=ax, cmap='viridis')
# Format
ax.set_axis_off()
Let's make it interactive!¶
Make sure to check the crs!
In [80]:
parks.crs
Out[80]:
The CRS is "EPSG:32618" and since it's not in 4326, we'll need to specify the "crs=" keyword when plotting!
In [81]:
# trim to only the columns we want to plot
cols = ["median_NDVI", "SITE_NAME", "geometry"]
# Plot the parks colored by median NDVI
p = parks[cols].hvplot.polygons(c="median_NDVI",
geo=True,
crs=32618,
cmap='viridis',
hover_cols=["SITE_NAME"])
# Plot the city limit boundary
cl = city_limits.hvplot.polygons(
geo=True,
crs=32618,
alpha=0,
line_alpha=1,
line_color="black",
hover=False,
width=700,
height=600,
)
# combine!
cl * p
Out[81]:
Exercise: Measure the median NDVI for each of Philadelphia's neighborhoods¶
- A GeoJSON file of Philly's neighborhoods is available in the "data/" folder
- Once you measure the median value for each neighborhood, make the following charts:
- A histogram of the median values for all neighborhoods
- Bar graphs with the neighborhoods with the 20 highest and 20 lowest values
- A choropleth showing the median values across the city
You're welcome to use matplotlib/geopandas, but encouraged to explore the new hvplot options!
In [82]:
landsat.crs.data['init']
Out[82]:
In [83]:
# Load the neighborhoods
hoods = gpd.read_file("./data/zillow_neighborhoods.geojson")
# Convert to the landsat
hoods = hoods.to_crs(landsat.crs.data['init'])
In [84]:
# Calculate the zonal statistics
stats_by_hood = zonal_stats(hoods, NDVI, affine=landsat.transform, stats=['mean', 'median'])
In [85]:
stats_by_hood
Out[85]:
In [86]:
# Store the median value in the neighborhood data frame
hoods['median_NDVI'] = [s['median'] for s in stats_by_hood]
In [87]:
hoods.head()
Out[87]:
Make a histogram of median values¶
In [88]:
# Initialize
fig, ax = plt.subplots(figsize=(8,6))
# Plot a quick histogram
ax.hist(hoods['median_NDVI'], bins='auto')
ax.axvline(x=0, c='k', lw=2)
# Format
ax.set_xlabel("Median NDVI", fontsize=18)
ax.set_ylabel("Number of Neighborhoods", fontsize=18);
Plot a (static) choropleth map¶
In [89]:
# Initialize
fig, ax = plt.subplots(figsize=(10,10))
# Plot the city limits
city_limits.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=4)
# Plot the median NDVI
hoods.plot(column='median_NDVI', legend=True, ax=ax)
# Format
ax.set_axis_off()
Plot an (interactive) choropleth map¶
With hvplot:
In [90]:
# trim to only the columns we want to plot
cols = ["median_NDVI", "ZillowName", "geometry"]
# Plot the parks colored by median NDVI
hoods[cols].hvplot.polygons(c="median_NDVI",
geo=True,
crs=32618,
width=700,
height=600,
cmap='viridis',
hover_cols=["ZillowName"])
Out[90]:
With altair:
In [97]:
# trim to only the columns we want to plot
cols = ["median_NDVI", "ZillowName", "geometry"]
(
alt.Chart(hoods[cols].to_crs(epsg=4326))
.mark_geoshape()
.encode(
color=alt.Color(
"median_NDVI",
scale=alt.Scale(scheme="viridis"),
legend=alt.Legend(title="Median NDVI"),
),
tooltip=["ZillowName", "median_NDVI"],
)
)
Out[97]:
Plot a bar chart of the neighborhoods with the top 20 largest median values¶
In [99]:
# calculate the top 20
cols = ['ZillowName', 'median_NDVI']
top20 = hoods[cols].sort_values('median_NDVI', ascending=False).head(n=20)
top20.hvplot.bar(x='ZillowName', y='median_NDVI', rot=90)
Out[99]:
With altair:
In [104]:
(
alt.Chart(top20)
.mark_bar()
.encode(x=alt.X("ZillowName", sort="-y"), # Sort in descending order by y
y="median_NDVI",
tooltip=['ZillowName', 'median_NDVI'])
.properties(width=800)
)
Out[104]:
Plot a bar chart of the neighborhoods with the bottom 20 largest median values¶
In [106]:
cols = ['ZillowName', 'median_NDVI']
bottom20 = hoods[cols].sort_values('median_NDVI', ascending=True).tail(n=20)
bottom20.hvplot.bar(x='ZillowName', y='median_NDVI', rot=90)
Out[106]:
With altair:
In [109]:
(
alt.Chart(bottom20)
.mark_bar()
.encode(x=alt.X("ZillowName", sort="y"), # Sort in ascending order by y
y="median_NDVI",
tooltip=['ZillowName', 'median_NDVI'])
.properties(width=800)
)
Out[109]:
That's it!¶
- Next week: introduction to APIs and downloading data from the web in Python
- See you on Monday!