Natural Constraints Post Processing

import geolime as geo
from pyproj import CRS, Transformer
import numpy as np
import osmnx as ox
import json

geo.Project().set_crs(CRS('EPSG:20350'))

bm = geo.read_file("../data/block_model.geo")

dh = geo.read_file("../data/dh_pop_classif.geo")

bm_fe_map = geo.plot_2d(bm, 'Fe_kriged', 'mean')

bm_fe_map

Road Filtering

centroid = np.mean(bm.bounds(), axis=0)

transform = Transformer.from_crs(geo.Project().crs.to_epsg(), CRS("epsg:4326")) 

lat, lon = transform.transform(centroid[0], centroid[1])

road_gdf = ox.geometries_from_point(
    (lat, lon), 
    dist=np.max(np.diff(bm.bounds(), axis=0)) / 2, 
    tags={'highway': ['secondary', 'track']} 
)

road_buff = road_gdf.to_crs(CRS('EPSG:20350')).buffer(100)

road_buff = geo.GISObject('road_buffer', road_buff)

/home/runner/work/rocklea-report/rocklea-report/jupylime/geolime/schema/object_db.py:839: FutureWarning:

You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.

road_buff.plot(interactive_map=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

geo.select_points_inside_polygon(road_buff, None, bm, 'in_road')

bm.delete_cells('in_road')

geo.plot_2d(bm, 'Fe_kriged', 'mean')

River Filtering

river_gdf = ox.geometries_from_point(
    (lat, lon), 
    dist=np.max(np.diff(bm.bounds(), axis=0)) / 2, 
    tags={'water': True} 
)

river_buff = river_gdf.to_crs(CRS('EPSG:20350')).buffer(100)

river_buff = geo.GISObject('river_buffer', river_buff)

/home/runner/work/rocklea-report/rocklea-report/jupylime/geolime/schema/object_db.py:839: FutureWarning:

You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.

river_buff.plot(interactive_map=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

geo.select_points_inside_polygon(river_buff, None, bm, 'in_river')

bm.delete_cells('in_river')

geo.plot_2d(bm, 'Fe_kriged', 'mean')

Lonely Blocks Removal

bm.delete_cells('Y < 7473193')

river_line_gdf = ox.geometries_from_point(
    (lat, lon), 
    dist=np.max(np.diff(bm.bounds(), axis=0)) / 2, 
    tags={'waterway': True} 
)

river_line_buff = river_line_gdf.to_crs(CRS('EPSG:20350')).buffer(2000, single_sided=True)

river_line_buff = geo.GISObject('river_line_buffer', river_line_buff)

/home/runner/work/rocklea-report/rocklea-report/jupylime/geolime/schema/object_db.py:839: FutureWarning:

You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.

river_line_buff.plot(interactive_map=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

geo.select_points_inside_polygon(river_line_buff, None, bm, 'in_extended_river')

bm.delete_cells('in_extended_river')

Create Global Map with every objects

bm_fe_map = geo.plot_2d(bm, 'Fe_kriged', 'mean')

bm_fe_map

bm_fe_map.add_choroplethmapbox(
    geojson=json.loads(river_buff._geometry._data.to_crs("EPSG:4326").to_json()), 
    locations=river_buff._geometry._data.index, 
    marker_line_width=0, 
    z=river_buff._geometry._data.index, 
    colorscale='Pubu'
)
bm_fe_map.add_choroplethmapbox(
    geojson=json.loads(road_buff._geometry._data.to_crs("EPSG:4326").to_json()), 
    locations=road_buff._geometry._data.index, 
    marker_line_width=0, 
    z=road_buff._geometry._data.index, 
    colorscale='Brwnyl'
)
bm_fe_map.update_layout(mapbox_style='white-bg')

bm_fe_map.update_traces(showscale=False)

bm.keep_only_cells('Ore')

m = bm.plot_2d("Fe_kriged", "mean", interactive_map=True, zoom_start = 15)

/home/runner/work/rocklea-report/rocklea-report/jupylime/geolime/schema/object_db.py:839: FutureWarning:

You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.

from folium.plugins import SideBySideLayers
import folium

# layer_right = folium.TileLayer('https://a.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png')
layer_left = folium.TileLayer(
        tiles = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}',
        attr = 'Esri',
        name = 'Esri Map',
        overlay = False,
        control = True
       )
layer_right = folium.TileLayer(
        tiles = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        attr = 'Esri',
        name = 'Esri Satellite',
        overlay = False,
        control = True
       )

sbs = SideBySideLayers(layer_left=layer_right, layer_right=layer_left)

layer_left.add_to(m)
layer_right.add_to(m)
sbs.add_to(m)

<folium.plugins.side_by_side.SideBySideLayers at 0x7fe6211ea700>

from branca.element import Figure
fig = Figure(width=1636, height=829)
fig.add_child(m)

fig.save('test.html')

## Use cloud convert to convert html to png