Natural Constraints Post Processing#

import geolime as geo
from pyproj import CRS, Transformer
import numpy as np
import osmnx as ox
import json
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 4
      2 from pyproj import CRS, Transformer
      3 import numpy as np
----> 4 import osmnx as ox
      5 import json

File /opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/__init__.py:3
      1 """OSMnx init."""
----> 3 from ._api import *
      4 from ._version import __version__

File /opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/_api.py:12
     10 from .elevation import add_node_elevations_google
     11 from .elevation import add_node_elevations_raster
---> 12 from .features import features_from_address
     13 from .features import features_from_bbox
     14 from .features import features_from_place

File /opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/features.py:19
     17 import geopandas as gpd
     18 import pandas as pd
---> 19 from shapely.errors import GEOSException
     20 from shapely.errors import TopologicalError
     21 from shapely.geometry import LineString

ImportError: cannot import name 'GEOSException' from 'shapely.errors' (/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/shapely/errors.py)
geo.Project().set_crs(CRS('EPSG:20350'))
bm = geo.read_file("../data/block_model.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)
road_buff.plot(interactive_map=True)
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)
river_buff.plot(interactive_map=True)
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)
river_line_buff.plot(interactive_map=True)
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)