Visualize data on an interactive map

This notebook uses lonboard for interactive visualisation of data.

import geopandas as gpd
import numpy as np
import lonboard
from core.utils import used_keys
from lonboard.colormap import apply_continuous_cmap
import matplotlib as mpl
from mapclassify import classify
from sidecar import Sidecar

Define data path

chars_dir = "/data/uscuni-eurofab-overture/processed_data/chars/"

Define region

region = 65806
region = 66292 # prague

Buildings

Load building data and ensure the geometries are all valid Polygons.

buildings = gpd.read_parquet(f"{chars_dir}buildings_chars_{region}.parquet").to_crs(4326).reset_index()

buildings.geometry = buildings.make_valid()

buildings = buildings[buildings.geom_type.str.contains("Polygon")]

Create a lonboard layer

%%time
layer = lonboard.PolygonLayer.from_geopandas(buildings, opacity=.3)

CPU times: user 2.47 s, sys: 239 ms, total: 2.7 s
Wall time: 2.7 s

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='buildings')

Create a Map object

m = lonboard.Map(layer)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

buildings.columns

Index(['index', 'id', 'source', 'class', 'x', 'y', 'iid', 'geometry', 'ssbCCo',
       'ssbCor', 'ssbSqu', 'ssbCCM', 'ssbCCD', 'sdbAre', 'sdbPer', 'sdbCoA',
       'ssbERI', 'ssbElo', 'stbOri', 'mtbSWR', 'libNCo', 'ldbPWL', 'mibCou',
       'mibLen', 'mibAre', 'mibElo', 'mibERI', 'mibCCo', 'mibLAL', 'mibFR',
       'mibSCo', 'ltcBuA', 'mtbAli', 'mtbNDi', 'ltbIBD', 'stbCeA', 'nID',
       'stbSAl', 'nodeID'],
      dtype='object')

Specify a column and pass its values into a choropleth representation within the map.

# buildings.explore()

column = 'mtbSWR'
used_keys[column]

'shared walls ratio of buildings'

classifier = classify(buildings[column], 'quantiles', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_fill_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])

/home/krasen/morphometrics/.pixi/envs/default/lib/python3.12/site-packages/mapclassify/classifiers.py:1653: UserWarning: Not enough unique values in array to form 20 classes. Setting k to 7.
  self.bins = quantile(y, k=k)

buildings[column].describe().iloc[1:]

mean    0.084658
std     0.175582
min     0.000000
25%     0.000000
50%     0.000000
75%     0.062194
max     1.000000
Name: mtbSWR, dtype: float64

import pandas as pd
trid = region
data = pd.read_parquet(f'{chars_dir}primary_chars_{trid}.parquet')
data_source = pd.read_parquet(f"{chars_dir}buildings_chars_{trid}.parquet", columns=['source', 'class'])
data['source'] = data_source['source']

-4519                         NaN
-4518                         NaN
-4517                         NaN
-4516                         NaN
-4515                         NaN
                    ...          
 638176             OpenStreetMap
 638177             OpenStreetMap
 638178             OpenStreetMap
 638179             OpenStreetMap
 638180    Microsoft ML Buildings
Name: source, Length: 642700, dtype: object

buildings.loc[[638179, 638180], 'source']

638179             OpenStreetMap
638180    Microsoft ML Buildings
Name: source, dtype: object

Tessellation

Load tessellation data and ensure the geometries are all valid Polygons.

tess = gpd.read_parquet(f"{chars_dir}tessellations_chars_{region}.parquet").to_crs(4326)
tess['geometry'] = tess.geometry.make_valid()

tess = tess[tess.geom_type.str.contains("Polygon")]

Create a lonboard layer

%%time
layer = lonboard.SolidPolygonLayer.from_geopandas(tess, opacity=.2)

CPU times: user 4.3 s, sys: 328 ms, total: 4.63 s
Wall time: 4.62 s

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='tess')

Create a Map object

m = lonboard.Map(layer)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

tess.columns

Index(['enclosure_index', 'geometry', 'stcOri', 'sdcLAL', 'sdcAre', 'sscCCo',
       'sscERI', 'mtcWNe', 'mdcAre', 'ltcWRB', 'sicCAR', 'stcSAl', 'nID',
       'nodeID'],
      dtype='object')

column = 'mtcWNe'
used_keys[column]

'perimeter-weighted neighbours of ETC'

Specify a column and pass its values into a choropleth representation within the map.

classifier = classify(tess[column], 'quantiles', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_fill_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])

Enclosures

Load data and ensure the geometries are all valid Polygons.

enc = gpd.read_parquet(f"{chars_dir}enclosures_chars_{region}.parquet").to_crs(4326)

enc.geometry = enc.make_valid()

enc = enc[enc.geom_type.str.contains("Polygon")]

Create a lonboard layer

%%time
layer = lonboard.SolidPolygonLayer.from_geopandas(enc, opacity=.3)

CPU times: user 168 ms, sys: 31.9 ms, total: 200 ms
Wall time: 200 ms

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='enclosures')

Create a Map object

m = lonboard.Map(layer)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

enc.columns

Index(['eID', 'geometry', 'ldkAre', 'ldkPer', 'lskCCo', 'lskERI', 'lskCWA',
       'ltkOri', 'ltkWNB', 'likWCe', 'likWBB'],
      dtype='object')

column = 'ldkPer'
used_keys[column]

'perimeter of enclosure'

Specify a column and pass its values into a choropleth representation within the map.

classifier = classify(enc[column], 'quantiles', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_fill_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])

Streets

Load data and ensure the geometries are all valid Polygons.

streets = gpd.read_parquet(f"{chars_dir}streets_chars_{region}.parquet")

streets.geometry = streets.make_valid()

Create a lonboard layer

%%time
layer = lonboard.PathLayer.from_geopandas(streets.to_crs(4326), width_min_pixels=1)

CPU times: user 223 ms, sys: 12.3 ms, total: 236 ms
Wall time: 235 ms

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='streets')

Create a Map object

m = lonboard.Map(layer)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

assert np.allclose(streets['sdsLen'] , streets.geometry.length)

if 'mm_len' in streets.columns:
    assert np.allclose(streets['mm_len'] , streets.geometry.length)

streets.columns

Index(['geometry', '_status', 'mm_len', 'cdsbool', 'node_start', 'node_end',
       'sdsLen', 'sssLin', 'ldsMSL', 'sdsAre', 'ldsRea', 'ldsAre', 'sisBpM',
       'sdsSPW', 'sdsSPO', 'sdsSWD', 'nID'],
      dtype='object')

streets[['sdsLen', 'sssLin', 'ldsMSL', 'sdsAre', 'ldsRea', 'ldsAre',
       'sisBpM', 'sdsSPW', 'sdsSPO', 'sdsSWD']].describe()

	sdsLen	sssLin	ldsMSL	sdsAre	ldsRea	ldsAre	sisBpM	sdsSPW	sdsSPO	sdsSWD
count	77374.000000	77374.000000	77374.000000	5.267900e+04	77374.000000	6.464400e+04	52397.000000	77374.000000	77374.000000	57729.000000
mean	183.432278	0.959614	196.258464	2.242239e+04	136.501137	3.738741e+05	0.064995	32.483240	0.702718	4.339666
std	363.929170	0.106236	170.058588	4.710648e+04	119.552318	3.704621e+05	0.071708	13.008305	0.250935	1.994679
min	2.000000	0.000000	3.741196	4.152421e-02	0.000000	4.002875e+01	0.000262	0.000000	0.000000	0.000000
25%	58.894980	0.977297	103.906101	3.563385e+03	40.000000	1.548616e+05	0.027475	21.761670	0.500000	3.055577
50%	97.801325	0.998925	134.558532	8.524029e+03	122.000000	2.476303e+05	0.052702	30.189255	0.701923	4.544631
75%	170.871458	1.000000	214.998976	2.084422e+04	202.000000	4.465998e+05	0.085624	50.000000	1.000000	5.754629
max	11193.707415	1.000000	5524.818599	1.276722e+06	1788.000000	4.554284e+06	3.977393	50.000000	1.000000	12.135564

Specify a column and pass its values into a choropleth representation within the map.

column = 'sdsSPW'
used_keys[column]

'width of street profile'

streets[column] = streets[column].fillna(0)

classifier = classify(streets[column].astype(int), 'equalinterval', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])

Nodes

Load data and ensure the geometries are all valid Polygons.

nodes = gpd.read_parquet(f"{chars_dir}nodes_chars_{region}.parquet").to_crs(4326)

Create a lonboard layer

%%time
layer = lonboard.ScatterplotLayer.from_geopandas(nodes, radius_min_pixels=2)

CPU times: user 69.3 ms, sys: 4.83 ms, total: 74.2 ms
Wall time: 73.6 ms

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='nodes')

Create a Map object

m = lonboard.Map(layer, basemap_style=lonboard.basemap.CartoBasemap.Positron)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

nodes.columns

Index(['x', 'y', 'mtdDeg', 'lcdMes', 'linP3W', 'linP4W', 'linPDE', 'lcnClo',
       'lddNDe', 'linWID', 'ldsCDL', 'xcnSCl', 'mtdMDi', 'nodeID', 'geometry',
       'sddAre', 'midRea', 'midAre'],
      dtype='object')

Specify a column and pass its values into a choropleth representation within the map.

nodes[['mtdDeg', 'lcdMes', 'linP3W', 'linP4W', 'linPDE', 'lcnClo',
       'lddNDe', 'linWID', 'ldsCDL', 'xcnSCl', 'mtdMDi', 'nodeID', 'geometry',
       'sddAre', 'midRea', 'midAre']].describe()

	mtdDeg	lcdMes	linP3W	linP4W	linPDE	lcnClo	lddNDe	linWID	ldsCDL	xcnSCl	mtdMDi	nodeID	sddAre	midRea	midAre
count	61838.000000	61838.000000	61838.000000	61838.000000	61838.000000	6.183800e+04	61838.000000	61838.000000	61838.000000	61838.000000	61838.000000	61838.00000	4.654500e+04	61838.000000	4.908000e+04
mean	2.567418	0.105057	0.683189	0.130161	0.183928	1.004460e-06	0.005973	0.009984	319.454293	0.027892	183.852532	30918.50000	2.587954e+04	20.845176	6.426385e+04
std	1.005752	0.060368	0.116292	0.098865	0.115463	6.709700e-07	0.007673	0.005341	612.723878	0.053212	280.458199	17851.23731	4.584542e+04	23.542853	8.924184e+04
min	1.000000	-1.000000	0.000000	0.000000	0.000000	0.000000e+00	0.000396	0.000000	0.000000	0.000000	2.311779	0.00000	1.159084e+00	0.000000	1.221930e+01
25%	1.000000	0.061856	0.625000	0.058824	0.111111	4.621507e-07	0.003492	0.005882	73.066534	0.000000	74.482961	15459.25000	5.936987e+03	3.000000	1.793014e+04
50%	3.000000	0.098592	0.695652	0.111111	0.169492	8.682536e-07	0.005624	0.009957	181.984598	0.000000	109.536100	30918.50000	1.170070e+04	16.000000	3.585703e+04
75%	3.000000	0.142857	0.757576	0.181818	0.238095	1.420323e-06	0.007255	0.013344	364.331816	0.047619	170.624712	46377.75000	2.549037e+04	31.000000	7.125945e+04
max	6.000000	0.800000	1.000000	0.777778	1.000000	7.283451e-06	0.534588	0.119944	22877.000260	1.000000	11193.902019	61837.00000	1.170909e+06	583.000000	1.806340e+06

column = 'mtdMDi'
used_keys[column]

'mean distance to neighbouring nodes of street network'

classifier = classify(nodes[column], 'quantiles', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_fill_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])

Compare buildings across the pipeline

import geopandas as gpd

# read all data
all_gdf = gpd.read_parquet(f'/data/uscuni-eurofab/microsoft_buildings/ms_ce.pq')

# read all regions
regions_datadir = "/data/uscuni-eurofab/"
region_hulls = gpd.read_parquet(
        regions_datadir + "regions/" + "ms_ce_region_hulls.parquet"
    )
region_hulls.shape

(474, 1)

# subset the data
%%time
bidxs = all_gdf.sindex.query(region_hulls.loc[65806, 'convex_hull'], predicate='intersects')

CPU times: user 12.8 s, sys: 1.83 s, total: 14.6 s
Wall time: 14 s

plotting = all_gdf.iloc[bidxs]

%%time
plotting['geometry'] = plotting.simplify(1)

CPU times: user 1.1 s, sys: 34 ms, total: 1.14 s
Wall time: 1.14 s

/home/krasen/morphometrics/.pixi/envs/default/lib/python3.12/site-packages/geopandas/geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)

from lonboard import SolidPolygonLayer, Map
from lonboard.basemap import CartoBasemap
from lonboard.colormap import apply_categorical_cmap
from palettable.colorbrewer.qualitative import Set3_12
from core.cluster_validation import get_color

layer = SolidPolygonLayer.from_geopandas(
    gdf=plotting, opacity=0.15
)

/home/krasen/morphometrics/.pixi/envs/default/lib/python3.12/site-packages/lonboard/_geoarrow/ops/reproject.py:97: UserWarning: Input being reprojected to EPSG:4326 CRS
  warnings.warn("Input being reprojected to EPSG:4326 CRS")

m = Map(layer, basemap_style=CartoBasemap.Positron)
m

region_id = 65806
regions_buildings_dir = '/data/uscuni-eurofab/regions/buildings/'
buildings = gpd.read_parquet(regions_buildings_dir + f'buildings_{region_id}.pq')

layer = SolidPolygonLayer.from_geopandas(
    gdf=buildings, opacity=0.15
)
m = Map(layer, basemap_style=CartoBasemap.Positron)
m

/home/krasen/eurofab_morphometrics/.pixi/envs/default/lib/python3.12/site-packages/lonboard/_geoarrow/ops/reproject.py:97: UserWarning: Input being reprojected to EPSG:4326 CRS
  warnings.warn("Input being reprojected to EPSG:4326 CRS")

region_id = 65806
buildings_dir = '/data/uscuni-eurofab/processed_data/simplified_buildings/'
buildings = gpd.read_parquet(buildings_dir + f'buildings_{region_id}.parquet')

layer = SolidPolygonLayer.from_geopandas(
    gdf=buildings, opacity=0.15
)
m = Map(layer, basemap_style=CartoBasemap.Positron)
m

/home/krasen/eurofab_morphometrics/.pixi/envs/default/lib/python3.12/site-packages/lonboard/_geoarrow/ops/reproject.py:97: UserWarning: Input being reprojected to EPSG:4326 CRS
  warnings.warn("Input being reprojected to EPSG:4326 CRS")

region_id = 65806
buildings_dir = '/data/uscuni-eurofab/processed_data/buildings/'
buildings = gpd.read_parquet(buildings_dir + f'buildings_{region_id}.parquet')

# layer = SolidPolygonLayer.from_geopandas(
#     gdf=buildings, opacity=0.15
# )
# m = Map(layer, basemap_style=CartoBasemap.Positron)
# m

Primary chars

import pandas as pd

region = 65806
tess = gpd.read_parquet(f"{chars_dir}tessellations_chars_{region}.parquet", columns = ['enclosure_index', 'geometry']).to_crs(4326)
tess['geometry'] = tess.geometry.make_valid()
tess = tess[tess.geom_type.str.contains("Polygon")]

char_type = 'primary'


if char_type == 'weighted_lag':
    data = pd.read_parquet(f'/data/uscuni-eurofab/processed_data/chars/lag_chars_{region}_{kernel}_{spatial_lag}_{bandwith_type}.parquet')
    data = data.astype(np.float32)
    data.index = data.index.astype(np.int32)
elif char_type == 'unweighted_lag':
    data = pd.read_parquet(f'/data/uscuni-eurofab/processed_data/chars/lag_chars_{region}_unweighted_{spatial_lag}.parquet')
    data = data.astype(np.float32)
    data.index = data.index.astype(np.int32)
elif char_type == 'primary':
    data = pd.read_parquet(f'{chars_dir}primary_chars_{region}.parquet')
else:
    raise Exception('Datatype does not exist')

plotting = pd.merge(tess, data, left_index=True, right_index=True, how='inner')

plotting = plotting[plotting.index > -1]

%%time
layer = lonboard.SolidPolygonLayer.from_geopandas(plotting, opacity=.2)

CPU times: user 5.03 s, sys: 349 ms, total: 5.38 s
Wall time: 5.38 s

Create a Sidecar view (assumes JupyterLab) for more comfortable experience.

sc = Sidecar(title='tess')

Create a Map object

m = lonboard.Map(layer)

Display map within the sidecar plugin

with sc:
    display(m)

List avaialable columns

tess.columns

Index(['enclosure_index', 'geometry'], dtype='object')

column = 'sdsSPW'
used_keys[column]

'width of street profile'

Specify a column and pass its values into a choropleth representation within the map.

classifier = classify(plotting[column], 'quantiles', k=20)
normalizer = mpl.colors.Normalize(0, classifier.bins.shape[0])
vals = normalizer(classifier.yb)
layer.get_fill_color = apply_continuous_cmap(vals, mpl.colormaps['viridis'])