OpenLISEM maps creation, Windows

README.md

1. Please setup all the needed gdal and pcraster tools in your PATH
2. Make sure you follow the file structure:

-- soils -- DEM -- roads -- scripts -- ....

3. IMPORTANT: find your ldd value yourself and change all the environment variables!!!

kampala.bat

SET DEFAULT_PROJECTION=EPSG:21096
SET DIR_LANDUSE=landuse2016
SET DIR_SOILS=soils
SET GPKG_LAYER_CHANNELS=infrastructure_2003_fixed
SET GPKG_LAYER_ROADS=roads_osm20181012
SET GPKG_FILENAME=./lisem/kampala.gpkg
SET DEM_FILENAME_IN=DEM\kampala_demARCs1.tif

:: change these values!!!
SET LDD_VALUE=1e9
SET DISCHARGE_POINT=19887
rem SET DISCHARGE_POINT=19696
SET OUTPUT_RESOLUTION=20
SET RESAMPLE_ARG=r4
SET INITIAL_MOISTURE_1=0.7
SET INITIAL_MOISTURE_2=0.7
SET DEM_FILENAME_OUT=demarc%OUTPUT_RESOLUTION%m.map
SET DEM_MASK_FILENAME_OUT=dem%OUTPUT_RESOLUTION%m.map
SET MASK_FILENAME_OUT=mask%OUTPUT_RESOLUTION%m.map

rem GOTO:skip_demarc
rem GOTO:skip_mask
rem GOTO:skip_soils
rem GOTO:skip_channels
rem GOTO:skip_roads
:: Selecting a catchment is done in PCRaster, so we have to convert this file
:: to PCRaster format (from .tif to .map), type in (GDAL is integrated in PCRaster to facilitate conversion):
pcrcalc demarc.map = "%DEM_FILENAME_IN%"

:: We will resample this to a lower resolution, from 5m to 20m, for the modelling, this takes some time!
resample demarc.map %DEM_FILENAME_OUT% -%RESAMPLE_ARG%

:: We also have to force the coordinate system from “bottom to top” because of a bug in PCRaster:
mapattr -s %DEM_FILENAME_OUT% -P yb2t
:skip_demarc

:: To determine the boundaries of a watershed, we first create a flow network, using the steepest
:: slope. In PCRaster this is called a Local Drain Direction network. The command to create a LDD
:: network is:
pcrcalc --lddin ldd.map = lddcreate( %DEM_FILENAME_OUT%, %LDD_VALUE%, %LDD_VALUE%, %LDD_VALUE%, %LDD_VALUE% )


:: We will now find a location that is a likely outlet for the catchment under investigation. First make
:: a cumulative flow network. This map has the value ‘1’ at the tops mof hills where the flow starts,
:: and each time you go a cell downslope a value ‘1’ is added. In the valley at the end of a network,
:: the value is the sum of all cells upstream.
pcrcalc ups.map = accuflux(ldd.map, 1)

:: Check the discharge point in the end of your study area in ups.map and set it as a DISCHARGE_POINT value.
:: It means that X upstream cells are connected to this cell.Creating the watershed (or catchment)
:: for this point is done as follows:
pcrcalc ws.map = catchment(ldd.map, if(ups.map eq %DISCHARGE_POINT% then nominal(1)))

:: and to select the watershed only and fill the surrounding are with missing value:
pcrcalc ws1.map = if(ws.map eq 1 then scalar(1))

:: finally we cut the research area down to size, removing all missing values to create the mask for the database:
resample ws1.map %MASK_FILENAME_OUT% -C

:: and we cut the dem down to size:
resample --clone %MASK_FILENAME_OUT% %DEM_FILENAME_OUT% %DEM_MASK_FILENAME_OUT%

pcrcalc %DEM_MASK_FILENAME_OUT% *= %MASK_FILENAME_OUT%

rem Cleanup...
rem del ws1.map ws.map ldd.map ups.map
:skip_mask
rem GOTO:EOF


rem Fix projection on vegetation, build and bare land
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% %DIR_LANDUSE%\frBU2016.tif %DIR_LANDUSE%\frBU2016.temp.tif
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% %DIR_LANDUSE%\frSOIL2016.tif %DIR_LANDUSE%\frSOIL2016.temp.tif 
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% %DIR_LANDUSE%\frVEG2016.tif %DIR_LANDUSE%\frVEG2016.temp.tif 

rem Mask vegetation, built and bare land
call scripts\resample_mask.bat %DIR_LANDUSE%\frBU2016.temp.tif %MASK_FILENAME_OUT% build20m.map
call scripts\resample_mask.bat %DIR_LANDUSE%\frSOIL2016.temp.tif %MASK_FILENAME_OUT% bare20m.map
call scripts\resample_mask.bat %DIR_LANDUSE%\frVEG2016.temp.tif %MASK_FILENAME_OUT% veg20m.map

rem Cleanup...
del %DIR_LANDUSE%\frBU2016.temp.tif %DIR_LANDUSE%\frSOIL2016.temp.tif %DIR_LANDUSE%\frVEG2016.temp.tif 

rem Mask over soil layer 1
call scripts\resample_mask.bat %DIR_SOILS%\bulkdens_sl2r.tif %MASK_FILENAME_OUT% bulkdens_sl2.map
call scripts\resample_mask.bat %DIR_SOILS%\clay_sl2r.tif %MASK_FILENAME_OUT% clay_sl2.map
call scripts\resample_mask.bat %DIR_SOILS%\orgcarb_sl2r.tif %MASK_FILENAME_OUT% orgcarb_sl2.map
call scripts\resample_mask.bat %DIR_SOILS%\sand_sl2r.tif %MASK_FILENAME_OUT% sand_sl2.map

rem Mask over soil layer 2
call scripts\resample_mask.bat %DIR_SOILS%\bulkdens_sl5r.tif %MASK_FILENAME_OUT% bulkdens_sl5.map
call scripts\resample_mask.bat %DIR_SOILS%\clay_sl5r.tif %MASK_FILENAME_OUT% clay_sl5.map
call scripts\resample_mask.bat %DIR_SOILS%\orgcarb_sl5r.tif %MASK_FILENAME_OUT% orgcarb_sl5.map
call scripts\resample_mask.bat %DIR_SOILS%\sand_sl5r.tif %MASK_FILENAME_OUT% sand_sl5.map

rem Produce saxtonpedotransfer maps
pcrcalc -f .\scripts\saxtonpedotransfer.mod sl2 1 %INITIAL_MOISTURE_1%
pcrcalc -f .\scripts\saxtonpedotransfer.mod sl5 2 %INITIAL_MOISTURE_2%

:skip_soils
rem GOTO:EOF

:: NOTE: it's impossible to create PCRaster maps directly from gdal_rasterize because of this error:
:: $ gdal_rasterize -of PCRaster -a_srs EPSG:21096 -tr 5 5 -a width -l roads_osm20181012 ./kampala.gpkg roads.temp.map
:: ERROR 1: PCRaster driver: attempt to create dataset with an illegal data type (Float64); use either Byte, Int32 or Float32.

rem Prepare channel data
gdal_rasterize -tr 20 20 -a width -l %GPKG_LAYER_CHANNELS% %GPKG_FILENAME% chanwidth.temp1.tif
gdal_rasterize -tr 20 20 -a bottom -l %GPKG_LAYER_CHANNELS% %GPKG_FILENAME% chandepth.temp1.tif
gdal_rasterize -tr 20 20 -a depthr -l %GPKG_LAYER_CHANNELS% %GPKG_FILENAME% chandepthr.temp1.tif

gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% chanwidth.temp1.tif chanwidth.temp.tif 
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% chandepth.temp1.tif chandepth.temp.tif 
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% chandepthr.temp1.tif chandepthr.temp.tif 

pcrcalc chanwidth.temp.map = chanwidth.temp.tif 
pcrcalc chandepth.temp.map = chandepth.temp.tif 
pcrcalc chandepthr.temp.map = chandepthr.temp.tif 

:: fix the bug in pcraster
mapattr -s chanwidth.temp.map -P yb2t
mapattr -s chandepth.temp.map -P yb2t
mapattr -s chandepthr.temp.map -P yb2t

:: in case of negative width
pcrcalc chanwidth.temp.map = max(chanwidth.temp.map, 0)

:: apply the mask
resample --clone %MASK_FILENAME_OUT% chanwidth.temp.map chanwidth20m.map
resample --clone %MASK_FILENAME_OUT% chandepth.temp.map chandepth20m.map
resample --clone %MASK_FILENAME_OUT% chandepthr.temp.map chandepthr20m.map

:: we need also a channel mask that has the value of 1 for the channel and 0 for the rest of the area:
pcrcalc chanmask20m.map = if (chanwidth20m.map ne 0, scalar(1), 0)

:: remove temp files
rem del chanwidth.temp.tif chandepth.temp.tif chandepthr.temp.tif chanwidth.temp1.tif chandepth.temp1.tif chandepthr.temp1.tif chanwidth.temp.map chandepth.temp.map chandepthr.temp.map

:skip_channels
rem GOTO:EOF

:: TODO fix problem with roads network, they are created empty
rem Prepare roads data
rem gdal_rasterize -tr 20 20 -a width -l %GPKG_LAYER_ROADS% %GPKG_FILENAME% roads.temp1.tif

rem pcrcalc roads.temp.map = roads.temp.tif 

rem :: fix the bug in pcraster
rem mapattr -s roads.temp.map -P yb2t

rem :: set zeroes where no data is available
rem pcrcalc roads.temp.map = if(roads.temp.map ne 0, roads.temp.map, 0)

rem :: apply the mask
rem resample --clone %MASK_FILENAME_OUT% roads.temp.map roads20m.map


gdal_rasterize -tr 5 5 -a ID roads/roads.shp roads.temp1.tif
gdalwarp -s_srs EPSG:32636 -t_srs %DEFAULT_PROJECTION% roads.temp1.tif roads.temp.tif 

pcrcalc roads.temp.map = roads.temp.tif

mapattr -s roads.temp.map -P yb2t

pcrcalc roads.temp.map = if(roads.temp.map ne 0 then scalar(20) else 0)

resample --clone %MASK_FILENAME_OUT% roads.temp.map roads20m.map

rem Cleanup...
del roads.temp.tif roads.temp1.tif roads.temp.map
:skip_roads
rem GOTO:EOF


pcrcalc zero.map = %MASK_FILENAME_OUT% * 0

pcrcalc -f .\scripts\lisem_20m.mod %DISCHARGE_POINT%

rem GOTO:EOF

resample_mask.bat

SET INPUT_FILENAME=%1
SET MASK_FILENAME=%2
SET OUTPUT_FILENAME=%3

: fancy way in windows to get the filename only...
for /F %%i in ("%INPUT_FILENAME%") do @SET FILENAME=%%~ni

pcrcalc %FILENAME%.map = %INPUT_FILENAME%
mapattr -s "%FILENAME%.map" -P yb2t
resample --clone %MASK_FILENAME% "%FILENAME%.map" %OUTPUT_FILENAME%

DEL "%FILENAME%.map"

Hi,

I'm currently working on a personal project that involves openlisem and I have been running into problems generating proper maps from open-sourced data. I happened to stumble upon this guide and I'm really impressed! I was hoping to connect and get some advice on generating these maps.

Hope to hear from you!
Cheers
Zile