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

new version openheid2indicator

parent e9e2339a
"""
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 gemiddelde openheid voor beide jaren per 500m grid cell
3. bereken verschil in gemiddelde openheid als percentage en als absoluut getal in ha
Analyse in overleg met WM en HM dd 23 feb 2021
"""
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 avgopen_18
zonal_stats = rast.zonal_stats(grid, open18.read(1), stats=['mean', 'count'], affine=open18.affine, nodata=open18.nodata)
gem18_ha = pd.Series(data=[x['mean'] for x in zonal_stats], index=grid.index)
count18 = pd.Series(data=[x['count'] for x in zonal_stats], index=grid.index)
# calculate avgopen_19
zonal_stats = rast.zonal_stats(grid, open19.read(1), stats=['mean', 'count'], affine=open19.affine, nodata=open19.nodata)
gem19_ha = pd.Series(data=[x['mean'] for x in zonal_stats], index=grid.index)
count19 = pd.Series(data=[x['count'] for x in zonal_stats], index=grid.index)
# Calcualte trend1819
diff_perc = np.subtract(np.multiply(gem19_ha.divide(gem18_ha), 100), 100)
# Absolute difference
diff_ha = gem19_ha.subtract(gem18_ha)
# Append to grid shapefile and write to file
grid['avgopen18'] = gem18_ha
grid['avgopen19'] = gem19_ha
grid['count18'] = count18
grid['count19'] = count19
grid['diff_perc'] = diff_perc
grid['diff_ha'] = diff_ha
grid.to_file(r'w:\projects\Landschapsmonitor\Openheid\e_indicator\cbsgrid\openheid_1819_v20210223.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