Project 4: Road Finding¶

by Lorn

Hypothesis¶

It is possible to extract roads from satellite imagery with enough clarity to distinguish them from non-road features.

Data¶

Massachusetts Roads. A dataset of RGB images and rasterized mask of the road segments in that image.

https://www.kaggle.com/datasets/balraj98/massachusetts-roads-dataset

Microsoft Building Footprints. Building footprints in vector data format.

https://planetarycomputer.microsoft.com/dataset/ms-buildings

0.5m DEM of Massachusetts from MassGIS.

https://www.mass.gov/orgs/massgis-bureau-of-geographic-information

In [1]:
import rasterio
import numpy as np
import pandas as pd
import warnings
from pathlib import Path
import geopandas as gpd
from shapely.geometry import box
import matplotlib.pyplot as plt
import folium
import ipywidgets as widgets
from IPython.display import display
import joblib
import ipywidgets as widgets
from IPython.display import display
from scipy.ndimage import sobel, gaussian_filter
from sklearn.ensemble import RandomForestClassifier
from skimage.exposure import rescale_intensity
from skimage.filters import frangi
from skimage.morphology import closing, dilation, rectangle as footprint_rectangle

import warnings
warnings.filterwarnings("ignore")

View #1¶

The first view is a display of the raw data. We'll be working with the MassRoads dataset primarily and joining in additional data when needed. The images in MassRoads are distributed over a variety of small and large cities in the state. Mostly urban and suburban areas.

In [32]:
DATA_DIR = Path("roads")  

records = []

for tif in DATA_DIR.rglob("*.tiff"):
    with rasterio.open(tif) as src:
        if not src.crs:
            continue
        
        left, bottom, right, top = src.bounds
        geom = box(left, bottom, right, top)

        records.append({
            "path": str(tif),
            "geometry": geom
        })

gdf = gpd.GeoDataFrame(records, crs="EPSG:26986")  
gdf_ll = gdf.to_crs("EPSG:4326")


center = [42.25, -71.8]

m = folium.Map(location=center, zoom_start=9, tiles="cartodb positron")

for _, row in gdf_ll.iterrows():
    geojson = folium.GeoJson(row["geometry"])
    geojson.add_to(m)

m
Out[32]:
Make this Notebook Trusted to load map: File -> Trust Notebook