Commit 38ff1b99 authored by Roelofsen, Hans's avatar Roelofsen, Hans
Browse files

diverse new scripts etc

parent 3a6699ea
......@@ -24,6 +24,7 @@ def create_batch_script(year):
'REPLACE_BB_FC': bd[year]['bb']['vlak'],
'REPLACE_T10_GDB': bd[year]['top10']['gdb'],
'REPLACE_SCRATCH': r'c:\apps\temp_geodata\openheid21\scratch\1_1_{}'.format(year)}
except KeyError as e:
print('invalid year: {}'.format(e))
sys.exit(1)
......
"""
Create point vector shapefile with points equally spaced (100m) over Netherlands, restricted to positive areas in BBG
mask and fitting within spatial extent of Viewscape rasters.
To be used as viewpoints input for Viewscape.
Hans Roelofsen, 28-09-2021
"""
import os
import geopandas as gp
import numpy as np
from shapely import geometry
import rasterio as rio
# Read raster with Bestand Bodemgebruik 2015 valid data (>0) mask
mask = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\c_viewscape_input\coderingen_en_mask\BBGMask\BBGExtent_100m.tif')
mask_arr = np.where(mask.read(1) > 0, True, False)
# Set extent of output shapefile in meters in RD-New
xmin = 10000
ymin = 300000
xmax = 280000
ymax = 625000
# Spacing between points
spacing = 100
# How many points?
ncols = int(((xmax - xmin)/spacing))
nrows = int(((ymax - ymin)/spacing))
# Sequence of x-coordinates
x_coords = [x for x in range(int(xmin + (spacing/2)), int(xmax - (spacing/2)+spacing), spacing)]
# Sequence of y-coordinates
y_coords = [y for y in range(int(ymin + (spacing/2)), int(ymax - (spacing/2)+spacing), spacing)]
# two ncol*nrow matrices with X and Y coordinates. See: https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html
xx, yy = np.meshgrid(x_coords, y_coords)
# Zip x and y coordinates together with boolean inside/outside mask. Create shapely Point when mask is True
points = [geometry.Point(x, y) for x, y, m in zip(np.flipud(xx).flatten(),
np.flipud(yy).flatten(),
mask_arr.flatten())
if m]
# Write to file
gp.GeoDataFrame(data={'id': [i for i,_ in enumerate(points)]},
geometry=points).set_crs("EPSG:28992").to_file(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\c_viewscape_input\viewpoints\pnts_100m_v2021.shp')
"""
Create grid of 500m cells, mirroring design of https://www.cbs.nl/nl-nl/longread/diversen/2021/statistische-gegevens-per-vierkant-en-postcode-2020-2019-2018/2-technische-gegevens
Hans Roelofsen, WEnR,
"""
import sys
sys.path.append(r'c:/apps/proj_code/benb_utils')
import benb
import geopandas as gp
import numpy as np
# extent in meters
left = 10000
bottom = 300000
right = 280000
top = 625000
size = 500
ncols = int(((right - left)/size))
nrows = int(((top - bottom)/size))
sq = benb.gen_squares(x_ul=left, y_ul=top, nrow=nrows, ncol=ncols, size=500)
# Add CBS codering
sq['C28992R500'] = ['E{:04d}N{:04d}'.format(int(e), int(n)) for e, n in zip(np.divide(sq.geometry.bounds.minx, 100),
np.divide(sq.geometry.bounds.miny, 100))]
# Spation join to provincies
provs = gp.read_file(r'c:\Users\roelo008\OneDrive - WageningenUR\b_geodata\provincies_2018\poly\provincies_2018.shp')
foo = gp.sjoin(sq, provs.loc[:, ['provincien', 'geometry']], how='left')
# Drop grid cells zonder Provincie
foo.dropna(subset=['provincien'], inplace=True)
# Some hokken may be coupled to > 1 provincie. Undo this
if len(set(foo.ID)) < foo.shape[0]:
foo.drop_duplicates(subset=['ID'], inplace=True)
foo.drop(labels=['ID', 'index_right', 'size'], axis=1).to_file(r'c:\apps\temp_geodata\CBS_500m_grid.shp')
# sq.to_file(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\CBS500m_grids\versieWENR\CBS_500m_grid.shp')
\ No newline at end of file
call c:\Users\roelo008\Miniconda3\Scripts\activate.bat
call conda activate mrt
SET gdal_dir=C:/Program Files/QGIS 3.10/apps/Python37/Scripts
SET zip_dir=w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20210929_pd01012018\1_zip
SET shp_dir=w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20210929_pd01012018\2_shp
SET tif_dir=w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20210929_pd01012018\3_tiff
SET naambasis=horizon_pd01012018
SET timestamp=%date:~0, 4%%date:~5,2%%date:~8,2%t%time:~0,2%%time:~3,2%%time:~6,2%
REM Merge shapefiles into single shapefile
python "%gdal_dir%"\ogrmerge.py -o %shp_dir%/%naambasis%_merged.shp %zip_dir%/*.shp -single
REM Waarom nog naar raster omzetten, als de rapportage plaatsvindt op basis van de CBS 500m grids?
REM Kan ook direct daarnaartoe samenvatten
REM Rasterize
gdal_rasterize -co COMPRESS=LZW -a_nodata 0 -a_srs epsg:28992 -a T_VISAREA -l %naambasis%_merged -te 10000 300000 280000 625000 -tr 100 100 %shp_dir%\%naambasis%_merged.shp %tif_dir%\%naambasis%_merged_%timestamp%.tif
REM Calculate to hectares
python "%gdal_dir%"/gdal_calc.py -A %tif_dir%\%naambasis%_merged_%timestamp%.tif --co COMPRESS=LZW --outfile=%tif_dir%/%naambasis%_merged_hectare_%timestamp%.tif --calc="A/10000"
REM copy template layer file and rename
Copy w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\b_layers\Openheidsklassn.qml %tif_dir%
W:
cd %tif_dir%
rename Openheidsklassn.qml %naambasis%_merged_hectare_%timestamp%.qml
......@@ -18,33 +18,61 @@ Analyse in overleg met WM en HM dd 23 feb 2021
"""
import os
import numpy as np
import pandas as pd
import rasterio as rio
import rasterstats as rast
import geopandas as gp
def classify(x, arr, class_names):
"""
Classify x into an array
:param x:
:param arr: array with class boundaries
:param class_names: class names
:return: class name
"""
if np.isnan(x):
return np.nan
else:
return class_names[np.digitize(x, arr)]
# 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')
grid = gp.read_file(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\CBS500m_grids\versieWENR\CBS_500m_grid.shp')
open18 = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20210929_pd01012018\3_tiff\horizon_pd01012018_merged_hectare_20210929t104342.tif')
open19 = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20210928_pd01012019\3_tiff\horizon_pd01012019_merged_hectare_20210928t140125.tif')
# calculate avgopen_18
zonal_stats = rast.zonal_stats(grid, open18.read(1), stats=['mean', 'count'], affine=open18.affine, nodata=open18.nodata)
zonal_stats = rast.zonal_stats(grid, open18.read(1), stats=['mean', 'count'], affine=open18.transform, 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)
zonal_stats = rast.zonal_stats(grid, open19.read(1), stats=['mean', 'count'], affine=open19.transform, 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
# Calculate trend1819
diff_perc = np.subtract(np.multiply(gem19_ha.divide(gem18_ha), 100), 100)
# Absolute difference
diff_ha = gem19_ha.subtract(gem18_ha)
# Classify difference in hectare to classes
arr_vals = np.array([-100, -50, -25, 25, 50, 100])
classes = ['Veel geslotener (< -100 ha)',
'Aanzienlijk geslotener (-50 - -100 ha',
'Beperkt geslotener (-25 - -50 ha)',
'Geen of nauwelijks verandering (-25 - 25 ha)',
'Beperkt opener (25 - 50 ha)',
'Aanzienlijk opener (50 - 100 ha)',
'Veel opener (>= 100 ha)']
diff_class = diff_ha.apply(classify, arr=arr_vals, class_names=classes)
# Append to grid shapefile and write to file
grid['avgopen18'] = gem18_ha
grid['avgopen19'] = gem19_ha
......@@ -52,4 +80,7 @@ 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')
grid['diff_class'] = diff_class
grid['src_2018'] = os.path.basename(open18.name)
grid['src_2019'] = os.path.basename(open19.name)
grid.to_file(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\e_indicator\vergelijking1819\openheid_1819_v20210828.shp')
......@@ -25,7 +25,7 @@ REM Top10 Directorie
SET T10Dir=REPLACE_T10_GDB
REM MAsk dir
SET maskdir=W:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\BBGMask
SET maskdir=w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\c_viewscape_input\coderingen_en_mask\BBGMask
SET scratch_dir=c:\apps\temp_geodata\openheid21\scratch\1_1_REPLACE_YEAR
SET out_dir=c:\apps\temp_geodata\openheid21\output
......
WORKFLOW for openheidskaart
door Hans Roelofsen, WEnR, september 2021
1. bepaal de peildatum, meestal 1 januari van een bepaald jaar
2. definieer de brondata behorende bij de peildatum in .\brondata.py
3. run `create_batch_script.py` met bronjaar als command line argument
4. run het zojuist gecreerde batch script om de opgaande elementenkaarten te maken
5. kopieer de opgaande elementenkaart naar ..\vs\DataFiles
6. genereer ViewScape runfiles met ..\vs\gen_runfiles.py De opgaande-elementenkaart moet hierin opgenomen worden als
`landscapegridname`. Andere opties zijn optioneel.
7. specificeer de directory met runfiles in ..\vs\Run.bat
8. draai run.bat. Viewscape gaat nu in batch draaien op Azure
9. kopieer de viescape Shapefile outputs vanaf Azure storage naar een lokale directory
10. update de locatie van de Shapefiles en de output locaties in .\postprocessing\merge_rasterize.bat
11. draai .\postprocessing\merge_rasterize.bat
12. draai .\postprocessing\openheid2indictorv2.py
13. haal koffie
14. klaar
"""
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')
"""
Relief indicator naar CLO tekst en data
Hans Roelofsen, WEnR, september 2021
Input CBS 500m grid (non-WEnR versie, dus met onduidelijke IDs) met relief getallen (zie: w:\PROJECTS\Landschapsmonitor\cIndicatoren\Relief\a_geodata\Format_Reliëf.pdf)
Provincie Shapefile
Analyse:
Per provincie: totaal hectare AGR_OPP_AANPERC (Oppervlakte AAN-percelen in deze CBS 500*500 rastercel)
percentage AGR_OPP_AANPERC ten opzichte van provincie oppervlak
totaal hectare AGR_OPP_MUTATIES_M2 (Som van de totale veranderingen in hoogte. Hierbij wordt de AAN-oppervlakte van de verschillende klassen (stijgingen en dalingen) gesommeerd)
percentage AGR_OPP_MUTATIES_M2 ten opzichte van provincie oppervlak
verhouding AGR_OPP_MUTATIES_M2/AGR_OPP_AANPERC in percentage
totaal hectare NAT_OPP_NATUURLIJKTERREIN (Oppervlakte natuurlijk terrein in NatuurBeheerPlannen van de provincies in deze CBS 500*500 rastercel)
percentage NAT_OPP_NATUURLIJKTERREIN ten opzichte van provincie oppervlak
totaal hectare NAT_OPP_MUTATIES_M2 (Som van de totale veranderingen in hoogte. Hierbij wordt de natuurlijk terrein-oppervlakte van de verschillende klassen (stijgingen en dalingen) gesommeerd)
percentage NAT_OPP_MUTATIES_M2 ten opzichte van provincie oppervlak
verhouding NAT_OPP_MUTATIES_M2/NAT_OPP_NATUURLIJKTERREIN
"""
import pandas as pd
import geopandas as gp
import os
import datetime
relief_gdb = r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Relief\a_geodata\IndicatorRelief_30_augo_2021\IndicatorRelief_20210325.gdb'
relief_fc = 'NL_500x500'
provs_src = r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\provincies\incl_water\provincies_2018.shp'
relief = gp.read_file(relief_gdb, layer=relief_fc)
provs = gp.read_file(provs_src)
provs.index = provs.provincien
# Koppel elke 500m cell aan een Provincie
foo = gp.sjoin(relief, provs.loc[:, ['provincien', 'geometry']], how='left')
# Drop grid cells zonder Provincie
foo.dropna(subset=['provincien'], inplace=True)
# Some hokken may be coupled to > 1 provincie. Undo this
if len(set(foo.VIERKANT_ID)) < foo.shape[0]:
foo.drop_duplicates(subset=['VIERKANT_ID'], inplace=True)
# Pivot to summarize per province
piv = pd.pivot_table(data=foo, index='provincien', values=['AGR_OPP_AANPERC', 'AGR_OPP_MUTATIES_M2', 'NAT_OPP_NATUURLIJKTERREIN',
'NAT_OPP_MUTATIES_M2'], aggfunc='sum')
# convert sq meters to hectares
piv = piv.divide(10000)
piv.rename(columns={'AGR_OPP_MUTATIES_M2': 'AGR_OPP_MUTATIES_HA',
'NAT_OPP_MUTATIES_M2': 'NAT_OPP_MUTATIES_HA',
'AGR_OPP_AANPERC': 'AGR_OPP_AANPERC_HA',
'NAT_OPP_NATUURLIJKTERREIN': 'NAT_OPP_NATUURLIJKTERREIN_HA'}, inplace=True)
# Join provincie areas in hectare
piv = piv.join(other=provs.area_ha, how='left')
# Bereken oppervlaktes in precentage fracties mutaties
piv['AGR_OPP_AANPERC_PERC'] = piv.AGR_OPP_AANPERC_HA.divide(piv.area_ha).multiply(100)
piv['AGR_OPP_MUTATIES_PERC'] = piv.AGR_OPP_MUTATIES_HA.divide(piv.area_ha).multiply(100)
piv['NAT_OPP_NATUURLIJKTERREIN_PERC'] = piv.NAT_OPP_NATUURLIJKTERREIN_HA.divide(piv.area_ha).multiply(100)
piv['NAT_OPP_MUTATIES_PERC'] = piv.NAT_OPP_MUTATIES_HA.divide(piv.area_ha).multiply(100)
piv['NIET_NAT_OF_AGR_PERC'] = 100 - (piv.AGR_OPP_AANPERC_PERC + piv.NAT_OPP_NATUURLIJKTERREIN_PERC)
piv['NIET_NAT_OF_AGR_HA'] = piv.area_ha - (piv.AGR_OPP_AANPERC_HA + piv.NAT_OPP_NATUURLIJKTERREIN_HA)
# Write to output excel
out_dir = r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Relief\h_compendium_teksten'
excel_out = r'CompendiumSamenvatting.xlsx'
with pd.ExcelWriter(os.path.join(out_dir, excel_out), mode='w') as writer:
metadata = {'sourcedata': os.path.join(relief_gdb, relief_fc),
'provincie source': provs_src,
'made by': os.environ.get('USERNAME'),
'made when': datetime.datetime.now().strftime('%d-%b-%Y_%H:%M:%S'),
'made with': }
pd.DataFrame(data=metadata, index=[0]).T.to_excel(writer, sheet_name='Metadata', index=True)
stats_ha.to_excel(writer, sheet_name='Statistiek_hectare', index=True)
stats_perc.to_excel(writer, sheet_name='Statistiek_precentage', index=True)
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