DEM Alkmaar
In [8]:
Copied!
import geopandas as gpd
import matplotlib.pyplot as plt
import rasterio
from rasterio.windows import from_bounds
from shapely.geometry import box
from waterlagen import datastore
from waterlagen.ahn.download import get_ahn_rasters
from waterlagen.ahn.interpolate import interpolate_ahn_tiles
from waterlagen.bag.download import get_bag_features
from waterlagen.bag.rasterize import rasterize_bag
from waterlagen.logger import init_logger
from waterlagen.raster import create_vrt_file
import geopandas as gpd
import matplotlib.pyplot as plt
import rasterio
from rasterio.windows import from_bounds
from shapely.geometry import box
from waterlagen import datastore
from waterlagen.ahn.download import get_ahn_rasters
from waterlagen.ahn.interpolate import interpolate_ahn_tiles
from waterlagen.bag.download import get_bag_features
from waterlagen.bag.rasterize import rasterize_bag
from waterlagen.logger import init_logger
from waterlagen.raster import create_vrt_file
Aanpak¶
We gaan een DEM maken voor een klein stukje Alkmaar in de volgende stappen:
- Downloaden AHN
- Downloaden BAG
- Dichtinterpoleren AHN
- Verrasteren BAG op basis van het AHN uit stap 1
- Raster 3 en 4 samenvoegen in 1 VRT
Hieronder specificeren we de logger; die schrijft direct een logfile in de DataStore
In [9]:
Copied!
name = "dem_alkmaar"
logger = init_logger(name=name, log_file=datastore.data_dir / f"{name}.log")
bbox = [110960, 515600, 112560, 516760]
poly_mask = box(*bbox)
gpd.GeoSeries([poly_mask], crs=28992).explore()
name = "dem_alkmaar"
logger = init_logger(name=name, log_file=datastore.data_dir / f"{name}.log")
bbox = [110960, 515600, 112560, 516760]
poly_mask = box(*bbox)
gpd.GeoSeries([poly_mask], crs=28992).explore()
Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Downloaden AHN¶
In [10]:
Copied!
ahn_vrt = get_ahn_rasters(poly_mask=poly_mask)
ahn_vrt = get_ahn_rasters(poly_mask=poly_mask)
Downloaden BAG¶
In [11]:
Copied!
bag_gdf = get_bag_features(bbox, layer="pand", source="wfs")
bag_gdf = get_bag_features(bbox, layer="pand", source="wfs")
INFO waterlagen.bag.download: Start downloading BAG bag:pand using WFS, 1000 features per page (download)
Downloading page: 07
Dichtinterpoleren AHN¶
In [12]:
Copied!
ahn_filled_vrt = interpolate_ahn_tiles(ahn_vrt_file=ahn_vrt)
ahn_filled_vrt = interpolate_ahn_tiles(ahn_vrt_file=ahn_vrt)
INFO waterlagen.ahn.interpolate: start interpolating M_19BZ1 (1/1) INFO waterlagen.ahn.interpolate: writing d:\repositories\waterlagen\data\processed_data\ahn_filled\M_19BZ1.tif INFO waterlagen.raster: VRT file created d:\repositories\waterlagen\data\processed_data\ahn_filled\ahn_filled.vrt
Verrasteren vloerpeilen¶
In [13]:
Copied!
out_tif = datastore.processed_data_dir / "bag_pand.tif"
rasterize_bag(
dem_raster=ahn_filled_vrt,
bag_gdf=bag_gdf,
bag_pand_tif=out_tif,
buffer_step_m=1,
)
out_tif = datastore.processed_data_dir / "bag_pand.tif"
rasterize_bag(
dem_raster=ahn_filled_vrt,
bag_gdf=bag_gdf,
bag_pand_tif=out_tif,
buffer_step_m=1,
)
Samenvoegen¶
In [14]:
Copied!
dem_vrt = datastore.processed_data_dir / "dem.vrt"
dem_vrt = create_vrt_file(
vrt_file=dem_vrt,
directory=[
ahn_filled_vrt.parent,
datastore.processed_data_dir,
],
)
with rasterio.open(dem_vrt) as src:
window = from_bounds(*bbox, transform=src.transform)
data = src.read(1, window=window, masked=True)
transform = src.window_transform(window)
fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(
data,
extent=(bbox[0], bbox[2], bbox[1], bbox[3]),
origin="upper",
cmap="terrain",
)
plt.axis("off")
plt.show()
dem_vrt = datastore.processed_data_dir / "dem.vrt"
dem_vrt = create_vrt_file(
vrt_file=dem_vrt,
directory=[
ahn_filled_vrt.parent,
datastore.processed_data_dir,
],
)
with rasterio.open(dem_vrt) as src:
window = from_bounds(*bbox, transform=src.transform)
data = src.read(1, window=window, masked=True)
transform = src.window_transform(window)
fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(
data,
extent=(bbox[0], bbox[2], bbox[1], bbox[3]),
origin="upper",
cmap="terrain",
)
plt.axis("off")
plt.show()
INFO waterlagen.raster: VRT file created d:\repositories\waterlagen\data\processed_data\dem.vrt