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.features_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)

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.features_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)

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.features_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)

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)

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

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