In [55]:

    

%matplotlib inline

from __future__ import division

import os
import fiona as fio
from matplotlib import pyplot as plt
import numpy as np
import rasterio as rio
from rasterio.features import rasterize
from shapely.geometry import asShape
from shapely.geometry import shape
from scipy import ndimage
from skimage.morphology import square
from skimage.morphology import opening
from skimage.morphology import watershed

import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 15, 15


    

In [5]:

    

dtm_path = 'sample-data/DTM.tif'
dem_path = 'sample-data/DEM.tif'
mask_path = 'sample-data/thematic-mask.tif'
naip_path = 'sample-data/NAIP.tif'

# Keys are paths and vals are buffer distances
thematic_mask_paths = {
    'sample-data/tl_2014_24043_areawater.geojson': 0,
    'sample-data/tl_2014_24_prisecroads.geojson': 3,
    'sample-data/tl_2014_54037_areawater.geojson': 0,
    'sample-data/tl_2014_54037_linearwater.geojson': 1,
    'sample-data/tl_2014_54037_roads.geojson': 3
}

# Define some additional commonly used variables
with rasterio.open(dtm_path) as src:
    dtm_meta = src.meta.copy()
    dtm_shape = src.shape
    dtm_nodata = src.nodata
    dtm_dtype = src.meta['dtype']
    dtm_affine = src.affine
    dtm_bounds = src.bounds
    dtm_x_min, dtm_y_min, dtm_x_max, dtm_y_max = dtm_bounds


    

In [15]:

    

# Generate a Normalized DTM by subtracting the DEM from the DTM to give
# us absolute height above ground while removing all but the
# most dense vegetation
with rio.open(dtm_path) as dtm, rio.open(dem_path) as dem:
    ndtm = dtm.read(indexes=1) - dem.read(indexes=1)


    

In [20]:

    

# Load the intersecting NAIP area
with rio.open(naip_path) as src:
    col_min, row_min = [int(i) for i in ~src.affine * (dtm_x_min, dtm_y_max)]
    col_max, row_max = [int(i) for i in ~src.affine * (dtm_x_max, dtm_y_min)]
    window = ((row_min, row_max), (col_min, col_max))
    naip = src.read(
        window=window,
        boundless=True,
        masked=True
    )


    

In [54]:

    

with rio.open(mask_path) as src:
    thematic_mask = src.read(indexes=1).astype(np.bool_)

plt.title("Thematic Mask")
plt.imshow(thematic_mask)
plt.set_cmap('binary')
plt.show()

plt.title("nDTM")
plt.imshow(ndtm)
plt.set_cmap('YlGnBu_r')
plt.show()

with rasterio.open(dtm_path) as dtm:
    plt.title("DTM")
    plt.imshow(dtm.read(indexes=1))
    plt.show()

with rasterio.open(dem_path) as dem:
    plt.title("DEM")
    plt.imshow(dem.read(indexes=1))
    plt.show()

plt.title("NAIP")
plt.imshow(np.transpose(naip, [1,2,0]))
plt.show()


    

In [56]:

    

# Buildings we're interested in are at least 2.5 meters tall.
# Mask everything else.
# This is a point that could use some work.  Implementing actual
# segmentation to group objects like bushes next to houses would
# yield a better elevation mask.
markers = np.zeros_like(ndtm)
markers[ndtm < 2.5] = 1
markers[ndtm >= 2.5] = 2

# Immediately subtract 1 to match other masks
elevation_mask = ~(watershed(ndtm, markers) - 1).astype(np.bool_)
plt.title("Elevation Mask")
plt.imshow(elevation_mask)
plt.show()


    

In [62]:

    

# The NAIP we're working with doesn't have NIR, which would be
# extremely useful for eliminating vegetation, so we have to compute
# a greenness index to use for the same purpose.
# g / sum(r, g, b)
green_mask = naip[1] / (
    naip[0].astype(np.float32) + 
    naip[1].astype(np.float32) + 
    naip[2].astype(np.float32)
)
green_mask = np.nan_to_num(green_mask)
green_mask = green_mask >= 0.38
plt.title("Green Mask")
plt.imshow(green_mask)
plt.show()

# Merge elevation, thematic, and green masks together and add black
# NAIP pixels and fill small non-masked areas.  The buildings we care
# about are at least 25 square meters.
merged_mask = ((naip[0] == 0) & (naip[1] == 0) & (naip[2] == 0)) \
    + thematic_mask + elevation_mask + green_mask
m_segments, m_segment_labels = ndimage.label(~merged_mask)
m_segment_sizes = np.bincount(m_segments.ravel()) > 25
m_segment_sizes[0] = 0
final_mask = ~(m_segment_sizes[m_segments])

# Quick morphological cleanup
# This could also be improved with some additional morphological operations
# but that is currently handled in vector space after polygonizing.
# It would be nice to do a more agressive opening and then re-segmenting
# with the new kernals for a morphologically better product.
final_mask = opening(final_mask.astype(np.bool_), square(3)).astype(np.bool_)
plt.imshow(final_mask)
plt.show()


    

    




/Users/wursterk/github/Geoprocessing-Examples/venv/lib/python3.4/site-packages/IPython/kernel/__main__.py:8: RuntimeWarning: invalid value encountered in true_divide






    













    

In [63]:

    

# Normally this would be where the output data would be vectorized
# and additional morphological and topological operations would be
# performed but for some reason the rasterio and Fiona drivers are
# clobbering each other so the results written back to a raster.
# meta = {
#     'driver': 'GeoJSON',
#     'schema': {
#         'geometry': 'Polygon',
#         'properties': {}
#     },
#     'crs': dtm_meta['crs']
# }
# with fio.open('Extracted-Buildings.geojson', 'w', **meta) as dst:
#     for geom, dn in rasterio.features.shapes(final_mask.astype(np.uint8), transform=dtm_affine):
#         dst.write(
#             {
#                 'type': 'Feature',
#                 'geometry': geom,
#                 'properties': {}
#             }
#         )
processed = []
for idx, (geom, dn) in enumerate(
    rasterio.features.shapes(final_mask.astype(np.uint8), transform=dtm_affine)):
    if dn == 0:
        geom = shape(geom)
        geom = geom.simplify(3)
#         geom = geom.buffer(1)  # Geometry is causing a segfault
#         geom = geom.buffer(-1)
        processed.append(geom)
rasterized = rasterize(
    shapes=processed,
    out_shape=ndtm.shape,
    transform=dtm_affine,
)
meta = dtm_meta.copy()
meta.update(
    COMPRESS='LZW',
    TILED='YES',
    dtype=np.uint8,
    nodata=None,
    count=1,
    affine=dtm_affine,
    transform=rio.guard_transform(dtm_affine)
)
with rasterio.open('Extracted-Buildings.tif', 'w', **meta) as dst:
    dst.write_band(1, rasterized.astype(dst.meta['dtype']))


    

In [ ]: