Commit e9e2339a authored by Roelofsen, Hans's avatar Roelofsen, Hans
Browse files

openheid indicator script

parent 80c1a58d
"""
Postprocessing of LGN 2018-2019 raster file.
Hans Roelofsen, 21-12-2020
"""
import geopandas as gp
import pandas as pd
import numpy as np
df = gp.read_file(r'\\wur\dfs-root\PROJECTS\Landschapsmonitor\Landgebruik\OneDrive_1_12-10-2020\LGN2018-LGN2019_5m_eliminated\bewerkt_HR\LGN2018-LGN2019_5m_eliminated.tif.vat.dbf')
lgn = pd.read_csv(r'\\wur\dfs-root\PROJECTS\Landschapsmonitor\Landgebruik\OneDrive_1_12-10-2020\LGN2018-LGN2019_5m_eliminated\bewerkt_HR\lgn_coderingen.csv')
code2desc = dict(zip(lgn.lgn_code, lgn.lgn_desc))
lgn_code2monitor_code = dict(zip(lgn.lgn_code, lgn.monitor_code))
lgn_code2monitor_desc = dict(zip(lgn.lgn_code, lgn.monitor_desc))
lgn_code2val_code = dict(zip(lgn.lgn_code, lgn.val_code))
df['lgndesc18'] = df.lgn18.map(code2desc)
df['lgndesc19'] = df.lgn19.map(code2desc)
df['moncode18'] = df.lgn18.map(lgn_code2monitor_code)
df['moncode19'] = df.lgn19.map(lgn_code2monitor_code)
df['mondesc18'] = df.lgn18.map(lgn_code2monitor_desc)
df['mondesc19'] = df.lgn19.map(lgn_code2monitor_desc)
df['valcode18'] = df.lgn18.map(lgn_code2val_code)
df['valcode19'] = df.lgn19.map(lgn_code2val_code)
df['mon1819dif'] = np.where(df.moncode18 == df.moncode19, 0, 1)
df['val1819dif'] = np.where(df.valcode18 == df.valcode19, 0, 1)
"""
Script om de openheidsberekeningen voor 2018 en 2019 samen te combineren in het 500m grid polygoon bestand
Input
z:\Landschapsmonitor\Openheid\d_viewscape_output\20200709_v2018\merged_tiff\horizonVS2018_merged_ha.tif
z:\Landschapsmonitor\Openheid\d_viewscape_output\20200709_v2019\merged_tiff\horizonVS2019_merged_ha.tif
z:\Landschapsmonitor\CBS500m_grids\shp\Vierkanten_Basis.shp
Merk op dat we **niet** horizonVS2018_merged_ha_v2.tif en horizonVS2019_merged_ha_v2.tif gebruiken. Deze grids
zijn net anders opgelijnd en passen niet nietjes in het 500m grid.
Analyse stappen
1. Lees grid shapefile
2. bereken SUM openheid 2018 per gridcel as sumopen18
3. bereken SUM openheid 2019 per gridcel as sumopen19
4. bereken trend as sumopen_19/sumopen_18 as trend1819
5. bereken indicator 0.9 > trend1819 > 1.1 --> 1 else 0 as indic1819
6. write grid shapefile w. additional columns to file
"""
import numpy as np
import pandas as pd
import rasterio as rio
import rasterstats as rast
import geopandas as gp
# read data
grid = gp.read_file(r'W:\projects\Landschapsmonitor\CBS500m_grids\shp\Vierkanten_Basis.shp')
open18 = rio.open(r'W:\projects\Landschapsmonitor\Openheid\d_viewscape_output\20200709_v2018\merged_tiff\horizonVS2018_merged_ha.tif')
open19 = rio.open(r'W:\projects\Landschapsmonitor\Openheid\d_viewscape_output\20200709_v2019\merged_tiff\horizonVS2019_merged_ha.tif')
# calculate sumopen_18
zonal_stats = rast.zonal_stats(grid, open18.read(1), stats='sum', affine=open18.affine, nodata=open18.nodata)
sumopen18 = pd.Series(data=[x['sum'] for x in zonal_stats], index=grid.index)
# calculate sumopen_19
zonal_stats = rast.zonal_stats(grid, open19.read(1), stats='sum', affine=open19.affine, nodata=open18.nodata)
sumopen19 = pd.Series(data=[x['sum'] for x in zonal_stats], index=grid.index)
# Calcualte trend1819
trend1819 = sumopen19.divide(sumopen18)
# Calculate indicator
indic1819 = pd.Series(np.where((trend1819 > 0.9) & (trend1819 < 1.1), 0, 1), index=grid.index)
# Append to grid shapefile and write to file
grid['sumopen18'] = sumopen18
grid['sumopen19'] = sumopen19
grid['trend1819'] = trend1819
grid['indic1819'] = indic1819
grid.to_file(r'w:\projects\Landschapsmonitor\Openheid\e_indicator\cbsgrid\openheid_20182019.shp')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment