Natural Constraints Post Processing

Natural Constraints Post Processing#

import geolime as geo
from pyproj import CRS, Transformer
import numpy as np
import osmnx as ox
import json
from folium.plugins import SideBySideLayers
import folium
from branca.element import Figure

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

bm_smu = geo.read_file("../data/block_model_smu.geo")

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

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

bm_fe_map

Road Filtering#

centroid = np.mean(bm_smu.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_smu.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:851: 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_smu, 'in_road')

bm_smu.delete_cells('in_road')

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

River Filtering#

river_gdf = ox.geometries_from_point(
    (lat, lon), 
    dist=np.max(np.diff(bm_smu.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:851: 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_smu, 'in_river')

bm_smu.delete_cells('in_river')

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

Lonely Blocks Removal#

bm_smu.delete_cells('Y < 7473193')

river_line_gdf = ox.geometries_from_point(
    (lat, lon), 
    dist=np.max(np.diff(bm_smu.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:851: 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_smu, 'in_extended_river')

bm_smu.delete_cells('in_extended_river')

Create Global Map with every objects#

bm_fe_map = geo.plot_2d(bm_smu, '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_smu.keep_only_cells('channel_hg')

m = bm_smu.plot_2d("Fe_kriged", "mean", interactive_map=True, zoom_start = 14)

/home/runner/work/rocklea-report/rocklea-report/jupylime/geolime/schema/object_db.py:851: 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.

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 0x7f98f42a7dd0>

fig = Figure(width=800, height=700)
fig.add_child(m)
Contents