Align aerial imagery from O-D points (current use case is mines and industry headquarters.

pair_locator.py

import geopandas as gpd
import pandas as pd
import contextily as ctx
from rasterio import open as rioopen
from rasterio import band
from rasterio.mask import mask
from rasterio.plot import reshape_as_image, show
from rasterio.transform import Affine
from rasterio.warp import calculate_default_transform, reproject, Resampling
import rioxarray
from numpy.ma import masked_where
from numpy import isnan, uint8, array, int16, amin, amax, arange, logical_or
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter
from shapely.geometry import box
import os
import requests
import utm

def wgs_to_utm(lat, lon):
    utm_z = str(utm.from_latlon(lat, lon)[2])
    if len(utm_z) == 1:
        utm_z = '0' + utm_z
    if lat >= 0:
        epsg = '326' + utm_z
    else:
        epsg = '327' + utm_z
    return epsg

def bounds(row, col, crs, radius):
    y, x = row[col].y, row[col].x
    epsg = wgs_to_utm(y, x)
    buff = gpd.GeoSeries([row[col]]).set_crs(epsg=4326).to_crs(epsg=epsg).buffer(radius).to_crs(epsg=crs)
    bounds = buff.total_bounds
    nsew = {
        'n': bounds[3],
        's': bounds[1],
        'e': bounds[2],
        'w': bounds[0]
    }
    return nsew

def project_raster(file, crs):
    dst_crs = f'EPSG:{crs}'
    with rioopen(file) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        kwargs = src.meta.copy()
        kwargs.update(
            crs = dst_crs,
            transform = transform,
            width = width,
            height = height
        )
        with rioopen(file, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=band(src, i),
                    destination=band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest
                )
    return None

def clip_n_scoot(file, buff, crs, base, radius, rnd=False):
    buff_dst_crs = buff.to_crs(epsg=crs)
    bounds_dst_crs = buff_dst_crs.total_bounds
    if rnd:
        shapes = [f for f in buff_dst_crs]            
    else:
        shapes = [box(*bounds_dst_crs)]
    with rioopen(file) as src:
        img, transform = mask(src, shapes, crop=True)
        
        if base is not None:
            with rioopen(base + '.tif') as base_src:
                meta = base_src.meta
                base_tl = base_src.transform * (0, 0)
                base_br = base_src.transform * (base_src.width, base_src.height)
                base_radius = abs(base_tl[1] - base_br[1]) / 2
                ratio = radius/base_radius
                scale = Affine.scale(ratio)
                translate = Affine.translation(
                    (base_src.width / 2), 
                    (base_src.height / 2)
                )
                meta.update(
                    count = src.count,
                    transform = base_src.transform * translate * scale * ~translate
                )
                if src.count > 1:
                    pass
                else:
                    meta.update(
                        dtype = int16
                    )
        else:
            meta = src.meta
            meta.update(
                driver = "GTiff",
                height = img.shape[1],
                width = img.shape[2],
                transform = transform
            )
    with rioopen(file, 'w+', **meta) as dst:
        dst.write(img)
    

def locate_that(row, column, name, radius = 4000, base = None, rnd = False):
    epsg = wgs_to_utm(row[column].y, row[column].x)
    buff = gpd.GeoSeries([row[column]]).set_crs(epsg=4326).to_crs(epsg=epsg).buffer(radius).to_crs(epsg=3857)
    bounds = buff.total_bounds
    
    file = name + '.tif'
    _ = ctx.bounds2raster(*bounds, file, source=ctx.providers.Esri.WorldImagery, zoom = 16)
    
    project_raster(file, epsg)
    clip_n_scoot(file, buff=buff, crs=epsg, rnd=rnd, base=base, radius=radius)
    
    return None

def align_grids(over_name, base_name):
    print("Aligning grids.")
    b = rioxarray.open_rasterio(base_name + '.tif')
    b.rio.to_raster(base_name + '.tif')
    o = rioxarray.open_rasterio(over_name + '.tif')
    o_matched = o.rio.reproject_match(b)
    o_matched.rio.to_raster(over_name + '.tif')
    return None

def get_dem(n, s, e, w, file = 'dem', upsample_factor = 2):
    endpoint = 'https://portal.opentopography.org/API/globaldem'
    q = {'demtype': 'SRTMGL1', 'north': n, 'south': s, 'east': e, 'west': w}
    try:
        print("Querying OpenTopography.")
        response = requests.get(endpoint, params = q)
        response.raise_for_status()
        with open(file + '.tif', 'wb') as dst:
            dst.write(response.content)
    except requests.exceptions.HTTPError as errh:
        print(errh)
    except requests.exceptions.ConnectionError as errc:
        print(errc)
    except requests.exceptions.Timeout as errt:
        print(errt)
    except requests.exceptions.RequestException as err:
        print(err)
        
def convert_gray(img):
    gray = 0.07 * img[:,:,2] + 0.72 * img[:,:,1] + 0.21 * img[:,:,0]
    img_gray = gray.astype(uint8)
    return img_gray

def map_plot(file, dims=(24,24), base = 'base', over = 'over', dem = 'dem', contour_smooth = 50):
    print("Plotting map.")
    fig, ax = plt.subplots(nrows = 1,  ncols = 1, figsize=(24,24))

    base_color = LinearSegmentedColormap.from_list('magenta_to_white', ['#ff0160', '#fff'])
    over_color = LinearSegmentedColormap.from_list('blue_to_white', ['#08306b', '#fff'])

    with rioopen(base + '.tif') as base_src:
        base = convert_gray(reshape_as_image(base_src.read()))
        ax.set_axis_off()
        show(
            base,
            ax=ax,
            cmap = base_color,
            transform = base_src.transform
        )

    with rioopen(over + '.tif') as over_src:
        over = over_src.read()
        over_gray = convert_gray(reshape_as_image(over))
        show(
            masked_where(
                logical_or(
                    over_gray==0, over_gray==255
                ), over_gray
            ), 
            ax=ax,
            cmap = over_color,
            transform = over_src.transform
        )
    
    with rioopen(dem + '.tif') as dem_src:
        dem = dem_src.read()
        masked_dem = masked_where(dem == -32768, dem)
        ndmin = amin(masked_dem)
        ndmax = amax(masked_dem)
        r = arange(ndmin, ndmax, 5)
        filtered = gaussian_filter(masked_where(masked_dem < ndmin, masked_dem), contour_smooth)
        ct = show(
            masked_where(
                filtered <= ndmin, 
                filtered
            ),
            ax=ax, 
            transform = dem_src.transform, 
            contour=True, 
            levels=r,
            colors='#fff',
            linewidths = 0.8,
            alpha = 0.6
        )
    fig.savefig(f'{file}.pdf', bbox_inches='tight', dpi=300)
    
def row_process(row, base_geom, base_size, base_id, over_geom, over_size, over_id, rnd=True):
    for tif in [f for f in os.listdir('./') if f.endswith('.tif')]:
        os.remove(os.path.join('./', tif))
    print("Fetching and processing base imagery.")
    base = locate_that(row, 
                       column = base_geom, 
                       radius = base_size, 
                       name = 'base'
                      )
    print("Fetching and processing overlay.")
    over = locate_that(row, 
                       column = over_geom, 
                       radius = over_size, 
                       base = 'base',
                       rnd = rnd, 
                       name = 'over'
                      )
    align_grids('over', 'base')
    print("Fetching and processing DEM.")
    # Download raster, write to dem.tif
    nsew = bounds(row, crs=4326, col=over_geom, radius=base_size + 500)
    get_dem(nsew['n'], nsew['s'], nsew['e'], nsew['w'], file='dem')
    # Project it to appropriate UTM zone.
    epsg = wgs_to_utm(array(nsew['n'], nsew['s']), array([nsew['w'], nsew['e']]))
    project_raster('dem.tif', epsg)
    # Create new buffer for clipping contours.
    buff = gpd.GeoSeries([row[over_geom]]).set_crs(epsg=4326).to_crs(epsg=epsg).buffer(base_size - ((base_size-over_size) / 2))
    clip_n_scoot('dem.tif', buff=buff, crs=epsg, rnd=True, base='base', radius=base_size - ((base_size-over_size) / 2))
    align_grids('dem', 'base')
    map_plot(f'{int(row[base_id])}-{int(row[over_id])}')
    return None
 
hq_df = gpd.read_file('./data/firms/hq_geography.geojson')
prop_df = gpd.read_file('./data/properties/geography.geojson')
own_df = pd.read_csv('./data/properties/ownership.csv')
overview_df = pd.read_csv('./data/properties/overview.csv')

df = hq_df.merge(own_df, left_on='id', right_on='CompanyID')
df = df.merge(prop_df, left_on='PropertyID', right_on='id')
df = df.merge(overview_df, left_on='PropertyID', right_on='PropertyID')
df = df.rename(columns={"geometry_x": "geometry_hq", "geometry_y": "geometry_prop"})
df = df[
    ((df["Asset Type"] == 'Mine') | (df["Asset Type"] == 'Mine Complex'))
    & (df["Work Type"] != 'Underground')
    & ((df["Development Status"] == 'Production') | (df["Asset Type"] == 'Closed') | (df["Asset Type"] == 'Decommissioning'))
].reset_index()

outputs = df.loc[8:9,:].apply(
    row_process,
    axis = 1, 
    base_geom = 'geometry_hq', 
    base_size = 2500, 
    base_id = 'CompanyID',
    over_geom = 'geometry_prop', 
    over_size = 1500,
    over_id = 'PropertyID'
)