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")
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:819: 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}
)
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/utils_geo.py:335: ShapelyDeprecationWarning:
Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/geometries.py:816: ShapelyDeprecationWarning:
__len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry.
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/geometries.py:816: ShapelyDeprecationWarning:
Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/utils_geo.py:426: ShapelyDeprecationWarning:
Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
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:819: 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}
)
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/utils_geo.py:335: ShapelyDeprecationWarning:
Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/osmnx/utils_geo.py:426: ShapelyDeprecationWarning:
Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
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:819: 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)