Geodata for palmpy

Below, it will be outlined how the corresponding vector and raster data shall look like to be able to be handled by the script. For each type of data a rough how-to procedure is given on how to perform certain operations in QGIS to reach the required state. It is assumed that a user with minor QGIS experience is able to follow the instructions perfectly fine. Note that the presented procedures are not the only way to achieve the desired outcome - other methods may include leveraging the powerful functions of the geopandas python library or any other GIS Software. It shall be ensured though that the end results fulfills the presented requirements for a usage in the script.

Quick Reference Table

The following table lists all attributes to the individual shapefiles that are needed, and which information needs to be included in the raster data. Entries in brackets are optional and as such not required to be added in order for palmpy to work. For more details to each file, click on the file topic and see, which namelist entries are required.

Feature Type	Format	Attribute 1	Attribute 2	Attribute 3	Attribute 4	Attribute 5
Satellite Image	Raster	RGB in 3 channels (standard, no mods needed)
Terrain DTM	Raster	Terrain height (zt) per pixel.
Land Cover (Bodenbedeckung, bb)	Polygon	OBJEKTART or value of src_pav_type in namelist
Paved Surfaces	Polygon	BELAGSART or value of src_pav_type in namelist
Building Footprints	Polygon	HEIGHT_TOP	HEIGHT_BOT	ID	BLDGTYP
Rastered 3d Vegetation	Raster	Top Height of tree in pixel
Single Trees	Polygon	HEIGHT_TOP	HEIGHT_BOT	ID	(LAI)	(RADIUS)
Tree Lines	Polygon	HEIGHT_TOP	HEIGHT_BOT	ID	(LAI)	(RADIUS)
Resolved Forests	Polygon	HEIGHT_TOP	HEIGHT_BOT	ID	(LAI)	(OBJEKTART)
Streets	Polygon	OBJEKTART	(RADIUS)	(STRASSENROUTE)	(STUFE)	(BELAGSART)
Crop fields	Polygon	OBJEKTART	(HEIGHT_TOP)	(CROP)	(ID)

*Note: 3D Vegetation can either be supplied as shp (Tree lines, Resolved Forests or Single Trees) or tiff (Rastered 3D vegetation).

Information: PALM requires raster data in the static driver file. The script will convert any non-raster data to raster format. The reason for working with vector data for preprocessing is that it is much simpler to edit a vector file than raster data.

Geodata Shopping List

This section outlines which data must be gathered in order to be able to perform the steps outlined in this documentation and to use the present code to its full potential. This will be presented on a per-variable basis.

Variable	Required Information and Format	Example Datasets
zt	DTM in raster format (geotiff, ascii) and sufficient resolution	swissALTI3D, NASA SRTM
vegetation_type	land cover dataset including natural surface types, vector format	swissTLM3D
water_type	land cover dataset including water surface types, vector format	swissTLM3D
pavement_type	dataset including pavement areas, or individual files for streets, railway lines, train stations, airfields, etc., with an indication of its surface type, vector format	swissTLM3D, OpenStreetMap
soil_type	information about soil types, vector format	?
surface_fraction	computed from vegetation, water and pavement datasets (either 0 or 1 currently possible in PALM)	-
lad	Single Trees, Tree Rows or Forest Information, including their height, width and shape/tree type. If these are missing, they can be possibly extracted from other datasets, vector (or raster) format	swissTLM3D, DTM, DSM
tree_id	Tree-ID is often included in the above mentioned datasets
buildings_2d	building footprints or 2.5D/LOD1 datasets. If not available, building heights could be extracted form other datasets., vector format	swissBUILDINGS3D, DSM, DTM, OpenStreetMap
buildings_3d	see buildings_2d, additionally the lower height of bridges, vector format	see above
vegetation_pars	Any further available information
building_pars	Any further available information
albedo_pars	Any further available information
orthoimage	Ortho-referenced satellite imagery, raster format	swissimage

If any of the above vector datasets is not available in vector format, the make_static script may be adapted to incorporate this fact (Hint: change rasterandcuttlm() to cutalti(), which lacks the rasterization step)

Geodata Sources for which palmpy is ready

You can check which source files you can readily process with palmpy by having a look into the dictfolder, which is in palmpy/staticcreation. In there the translation dictionaries are stored. Readily usable are the following (section is work in progress):

Geodata Source	Notes	Links
Datenmodell Kanton Zürich DM01AVZH24	This is open data available from the Amtlichen Vermessung of Kanton Zürich in Switzerland. Similar data might be available for other Cantons of Switzerland with minor differences. Check first if the "translation" into palm categories makes sense by checking the file in the dict folder.	City of Zurich
Downscaled Switzerland Land Use/Cover Data LULC_down_unige_2018	This public dataset contains a raster of switzerland in 25m resolution and maps tiles to land use classes (about 80?). Suitable for simulations of mountainous terrain. Contains also building pixels, but mapping dict does not treat those! Make sure you check if a lot of houses are in your domain, and if they are, treat them with a shapefile or with a suitable fillvalue. More guidance and testing to follow.	Paper
TLM	Used for swissTLM3D datasets. Rather coarse shapefile (opendata) provides information about Land use in Switzerland.	swisstopo
Tutorial	Used for tutorial case, but basically DM01AVZH24	-

A note on Coordinate Reference Systems (CRS)

TLDR: Transform all your geodata into a single conformal CRS before proceeding to create a static driver from it. Choose a CRS with an EPSG Code (e.g. for Switzerland, LV95 is EPSG code 2056, LV03 is 21781.

At this stage, an important note needs to be brought up regarding coordinate reference systems (CRS) and the impact they have on the whole static file generation process. While users with some familiarity with geodata will be comfortable with CRS, but to new users this could prove to be an obstacle in the learning process.

Geodata is basically data that is linked to certain geographical coordinates. While for raster layers this geographical data may be fixed by an origin and certain extents on the globe, it could be for vector layers that each point has its own geographical coordinate. These geographical coordinates are not of much use without knowing which geographical map projection they refer to. There are multiple map projections available, each with its own strengths and weaknesses. Some may be conformal (true angles) and others equal-area (true areas), and depending on the type of question that is asked, another CRS may be more appropriate.

PALM calculates all its equation on a cartesian grid, where the axes are perpendicular amongst each other. Therefore, it is recommended to transform any geodata into a CRS that has somewhat equal properties, should reality be reproduced correctly.

For Switzerland, source data from swisstopo are mostly available in either Switzerlands own CRS LV03 or LV95, which can be directly processed with the static creation script. Should locations out of Switzerland be of interest, geodata could be transformed into an EPSG coordinate system. As it is not possible to represent the whole world with a single, large conformal CRS, there is a regional EPSG coordinate system available for practically every location on earth. These EPSG CRS have a validity range, i.e. outside of these valid extents the errors get rather large, so another EPSG CRS may be more appropriate. Another possibility is the usage of UTM coordinates.

xkcd comic #2256

Rastered Data

The static generation script requires some data to be present in rastered form, i.e. as a geotiff. For the moment, the format needs to be geotiff (a .tif-File with a header containing geographical metadata) - other formats may be possible in the future. The following data can be supplied:

• Orthoimage
• Digital Surface/Height Model (Digital Terrain Model DTM (german: Digitales Höhenmodell DHM or Digitales Geländemodell DGM)) (Note: Topography only, without buildings or trees!)

optional:

• Resolved vegetation from a digital surface model ( Digital Surface Model DSM (german: Digitales Oberflächenmodell DOM)). This file should only contain trees with their height relative to the ground. This could come from LIDAR-scans.

Orthoimage

An orthoimage refers to a georeferenced image of the area in question. In the script, if an orthoimage is supplied, it is clipped to the domain extents and nest and proble locations will be plotted on it. For an example, refer to the following image.

Plot of selected domain extents and probe locations on the clipped orthoimage.

Topography Data

• relevant namelist flag: doterrain

Topography data is usually available in raster format. Notable examples are swisstopo's swissALTI3D DTM dataset with two meter resolution or NASA's SRTM Dataset with one arcsecond resolution (about 20x30m (N/E) in swiss latitudes). Note that there are various represenations of topography data available. To create a PALM static driver file, we require a DTM, a digital terrain model, which includes the height of the underlying topography only. Often, there is also a DSM, a digital surface model, available, which includes features such as trees and buildings in the dataset. It is assumed that building information comes from separate data files, hence a DTM is intended to be used in conjunction with the static generation script.

An extract of Switzerland from the NASA SRTM DTM dataset.

Optional: Resolved Vegetation

• relevant namelist flag: veg_format in settings section, veg_height_bot

From May 2022, it is possible to provide high resolution vegetation as raster into palmpy. For that, a digital surface model is needed that only contains tree height information. The tree crown bottom height is given by the parameter veg_height_bot. In the city of Zurich, this height is usually 3 meters. The path to the raster file is given in the paths section in the namelist, under veg_raster. The alpha, beta and LAI values for each individual domain are taken from the parameters lai_forest, a_forest, b_forest.

Example for rastered tree heights.

Vector Data - Shape Files

While raster data requirements are straightforward, it gets more complicated with data in vector format. This is due to the manifold ways information can be included in a shapefile. There are different types of features (line, points, polygons) and the naming of fields (columns in an attribute table) are more or less up to the entity that produces the data (unproven claim by the author based on practical experience with multiple datasets).

Note: A shape file may consist of up to six files with different file endings that separately contain information like attributes, coordinate reference systems, etc. When moving files around on the computer, move all files together. When providing files in a namelist, provide the filename with .shp in the end.

In order for the workflow to work, the vector data supplied to the script needs to conform to certain requirements, namely it shall possess certain fields with a certain title and data in a certain format. An overview over the required attribute fields for each file, refer to the following table.

Land Cover (Bodenbedeckung)

• relevant namelist flag: dolandcover

Land Cover in the Glarner Sernftal from the swissTLM3D BB dataset, overlay with a satellite image. Note that not all surfaces are covered!

• Required feature type: Polygon

• Required attributes: OBJEKTART

This file shall contain land cover information including vegetation and water surfaces. The focus here lies on natural land cover, such as rocks, desert, glaciers etc. Artificially paved or sealed areas are to be put in a separate file.

If features are added manually, make sure the classification matches the mapping dictionary or assign a new class number that is also included as a new entry in the mapping dictionaries. A good practice is to give manually defined classes an identifier starting with 100X or similar, so they can be clearly identified as individual changes to the dataset.

This file (and the crops dataset described further below) will decide, which pixels in the static driver will be classified as a vegetation_type and which ones as a water_type. Pixels that are not covered by a polygon in those datasets will be set to a bulk value that is to be set in the namelist. Paved or sealed surfaces will be added on top of this information and overwrite any set vegetation_type or water_type.

Paved or Sealed Surfaces

relevant namelist flag: dopavedbb

Basic paved land cover around Zurich airport.

• Required feature type: Polygon

• Required attributes: BELAGSART

• Optional/Planned attributes: STUFE, OBJEKTART, STRASSENROUTE

While there may be specific datasets available that contain the necessary information to construct the land cover regarding water surfaces (most obvious) or vegetation, there is hardly any dataset available that specifically contains paved and sealed surfaces (Note: Personal experience by the author). Based on the employed pavement surface classification in PALM, information about asphalt, concrete, sett, gravel, cobblestone and other surfaces of that kind need to be collected in this single file. As the source of this information is most definitely heterogeneous as well, specific guidance how to proceed cannot be given in this location. It is recommended to gather all information available that may relate to paved or sealed surfaces, convert them if necessary to a polygon shape file and merge all information into a single file.

For example, working with swissTLM3D data, to get an initial pavement/sealed surface dataset, one could collect streets, traffic areas ("Verkehrsareale") and railway lines (the underlying gravel respectively) information and summarize it in a single file. Processing all those single datasets can be done applied with the described operations in QGIS.

Buildings

Buildings footprints (orange) of the Zürich Zoo area, including the Masoala Halle, the Elephant Building in the north and the FIFA building in the south with various sport areas around it (green).

• Required feature type: Polygon
• Required attributes: ID, HEIGHT_TOP, HEIGHT_BOT, BLDGTYP
• Optional/Planned attributes:

In order to perform urban climate simulations, the positions and heights of buildings is an important metric to include in a simulation. Therefore, we need to create a shape file that contains these two pieces of information. This mode of representation has a downside as well. Buildings with round roof shapes need to represented by multiple polygons, if the roof shape shall be represented in a realistic way. Looking at the Masoala Halle at the Zurich Zoo with is rounded shape, if it is represented with a single polygon as in the image above, the Masoala Halle will only be represented by its mean height in a simulation. For coarse resolutions, this approach is definitely fine. However, if resolutions approach 4 or 2 meters, complex roof shapes can be resolved and may influence the flow also in a different way than a flat roof would. This must be kept in mind when creating building information.

The source information to create a building dataset may come in various forms. While for example the City of Zurich publicly shares a 2.5D shapefile, which contains all building footprints and building height information, for some other regions this will be more complicated. For example may information only be available as a dxf file or an stl file, from where new procedures for data preprocessing need to be found. Also there is a possibility for building information to be available in multipolygon format, which could be processed in ArcPro from ESRI, but not with open source tools.

Depending on the available source data, construct a building dataset containing information about its height, its lower limit in space, its type and an identification number. The lower limit is used to represent bridges or overhanging buildings in the dataset. The identification number will allow individual changes to buildings based on their identification numbers. A method for this is not implemented though.

Note: Check the source, where the building height is measured from and compare it with the PALM definition of building height. As an example, PALM bulidings are put on top of the highest topo grid point in the building footprint pixels. That is quite opposite to e.g. the definition for building height in the Zurich LOD1 building dataset, where a the height represents the height of the building measured from the lowest intersection with topography, plus three meters to account for any underground stories. A possible workaround could be to do a zonal statistics for each building polygon to get min and max of the topography for each building. Then subtract the spread from the building height field. Also subtract three meters for the underground stories (not relevant in PALM simulation, at least right now).

Possible workflow to forge a simple building footprint shapefile paired with LIDAR DTM-DSM data into a compliant buildings dataset:

1. get rid of obsolete data fields, as they may slow QGIS down significantly.
2. Create new fields HEIGHT_TOP, HEIGHT_BOT, BLDGTYP and ID, if not yet present.
3. Assign a building type based on the building ID. If not available, assign bulk values.
4. Perform zonal statistics to assign building heights to the features.
5. Use the calculated median value as HEIGHT_TOP attribute.
6. Check, if the performed classifications and rooftop heights are realistic.

Resolved Vegetation

Since May 2022, it is possible to supply resolved vegetation as a raster file. Refer to the 'Rastered Data' section above.

• relevant namelist flag: dolad

In PALM, it is possible to resolve vegetation instead of just parameterizing it. Resolving vegetation acts as a momentum sink for the flow, affects radiation calculations and influences latent heat balances. The script requires information about vegetation positions and heights and forges this information into a 3D leaf area density array. The vertical leaf area density profile is calculated based on chosen values for the leaf area index (LAI) and shape parameters alpha and beta. The LAI values are currently provided as bulk values for all single trees, all tree lines and all forests separately. It is planned to have the LAI as an attribute in the corresponding shape files to allow for a heterogeneous definition of this parameter to account for different types of trees.

The script know where to place vegetation based on polygons provided by one to three vegetation files. The script also collects the provided vegetation information and turns it into a leaf area density array based on the supplied field values for HEIGHT_TOP and HEIGHT_BOT and some shape parameters. Assuming the tree top view corresponds to the polygon information provided to the script, the tree is rastered and will be represented as follows:

Illustration how a tree polygon is rastered. View from above and from the side.

The leaf area density array will be filled between HEIGHT_BOT and HEIGHT_TOP based on the chosen LAI and two shape parameters alpha and beta. The following figure provides an overview about possible tree shapes that can be constructed with alpha and beta. This figure can be called with the function call mst.showbetadistribution() with a loaded palmpy module.

Resulting shapes for different alpha and beta values. A fractional tree height of 1 represents tree top, 0 the HEIGHT_BOT value, so the lower limit of the tree crown.

With that, it is possible to create a variety of tree crowns. Currently, only the leaf area density is used in a palm simulation. In the future, also the basal area density of a tree will be used in calculation. This feature is not yet implemented in the script.

There is the possibility to provide resolved vegetation information with three files: one for resolved forest patches, one for tree lines and one for single trees. This approach emerged from the fact that forest patches are usually present as polygons, tree lines as line information and single trees as point information. In order to be processed in the script, these all need to be transformed into polygon features, which requires different treatments of the different source features types. If desired, all information can be provided in a single file (leave the other path variables in the namelist blank).

Single Trees

Single Trees at the lakefront in Yverdon.

• Required feature type: Polygon

• Required attributes: HEIGHT_TOP, HEIGHT_BOT, ID

• Optional/Planned attributes: RADIUS, LAI

Single trees are mostly available as point features in a shape file. We need to forge this into a polygon shape file. To that end, the function "Buffer" (german: Puffern) can be used, which creates a circular polygon around the point feature with a defined radius. The radius can be defined for each point feature individually, which is why I recommend to include a RADIUS attribute for each point. The radius can be set to a bulk value and changed individually at known locations, or the information can be copied from another field if available.

The height of a single tree is represented by the fields HEIGHT_TOP and can also be set to a bulk value. Alternatively, if a LIDAR DTM and DSM is available, tree height can also be extracted by means of a zonal statistics. The parameter HEIGHT_BOT represents the bottom limit of the tree crown. Between those values the crown / leaf area density array is constructed corresponding the shape parameters alpha and beta.

Possible workflow coming from a point shape file:

1. create RADIUS attribute with bulk value
2. modify RADIUS where tree widths are either relevant, exceptionally large or small
3. Puffer by RADIUS: creates polygon shape file
4. Add new attributes HEIGHT_TOP and HEIGHT_BOT
5. 1. assign bulk values and be done
   2. perform Zonal statistics (requires LIDAR DTM minus DSM data (so trees heights are present))
6. Check if the performed actions make sense and are realistic.

Tree Lines

Tree lines and single trees around Yverdon airfield.

• Required feature type: Polygon

• Required attributes: HEIGHT_TOP, HEIGHT_BOT, ID

• Optional/Planned attributes: RADIUS, LAI

Tree lines are often available as line features in a shape file. We need to transform this into a polygon shape file. The same techniques as for single trees can be applied. Assign a RADIUS (if needed based on the OBJEKTART, which could point to bushes, hedges or actual tree lines) to each line and Buffer it to get a Polygon shape file. Add HEIGHT_TOP and HEIGHT_BOT attributes as well, populate them with either bulk values based on OBJEKTART or perform zonal statistics again.

Possible workflow from tree line shape file to a polygon shape file:

1. create RADIUS attribute with bulk value based on the OBJEKTART with the field calculator
2. modify RADIUS where tree widths are either relevant, exceptionally large or small
3. Puffer by RADIUS: creates polygon shape file
4. Add new attributes HEIGHT_TOP and HEIGHT_BOT
5. 1. assign bulk values and be done
   2. perform Zonal statistics (requires LIDAR DTM minus DSM data (so trees heights are present))
6. Check if the performed actions make sense and are realistic.

Resolved Forests

A resolved forest with some visible tree lines and single trees north of the Yverdon airfield.

• Required feature type: Polygon
• Required attributes: HEIGHT_TOP, HEIGHT_BOT, ID
• Optional/Planned attributes: OBJEKTART, LAI

Larger vegetation patches, such as forests or grass fields that are higher than one or multiple grid cells can be present as polygons in the source datasets. As they are already represented as polygons, there is no need for buffering or similar techniques. However, the information may be extracted from a dataset with more than just the required information. This can be done by copying the shape file and deleting the obsolete data, in order to only retain the relevant vegetation classes. For example, the swissTLM3D BB dataset contains three types of land cover classes that can be represented by resolved vegetation: Gebüschwald/Brush Land (OBJEKTART 6), Forest (OBJEKTART 12) and interrupted forest (OBJEKTART 13). This information can be subset into a new file. Again, HEIGHT_TOP and HEIGHT_BOT can be gotten by applying zonal statistics using a LIDAR DTM-DSM dataset.

Possible workflow from a swissTLM BB polygon shape file to a finished resolved forests shape file:

1. Select Objects by expression (Ctrl+F3)
2. Create query for OBJEKTART != 6 AND OBJEKTART != 12 AND OBJEKTART !=13. This selects all objects that are not vegetation.
3. Delete the selected features.
4. Add fields for HEIGHT_TOP and HEIGHT_BOT
5. 1. Use the Field calculator to assign different bulk values for the height fields based on OBJEKTART.
   2. perform Zonal statistics (requires LIDAR DTM minus DSM data (so trees heights are available))
6. Check if the performed actions make sense and are realistic.

Street Types (for Emission Modeling)

• relevant namelist flag: dostreettypes

Required streets dataset as polygons (extract south of Zurich airport).

• Required feature type: Polygon
• Required attributes: OBJEKTART
• Optional/Planned attributes: RADIUS, BELAGSART, STRASSENROUTE, STUFE

Streets appear in the static driver twofold: once as an area of paved or sealed surfaces and once as dedicated street type array, which classifies streets by their relevance for emission contributions. The latter is only required if PALM is run with chemistry activated. As street information is often available as line features, they need to be buffered again.

Practical note: For simulations in Switzerland, streets could be taken from Openstreetmap data or from swissTLM3D data. The OSM street dataset contains a clear classification by street type (highway, trunk, residential etc.), which the PALM street classification system is also built upon. However, The OSM dataset does not include a classification by width. The swissTLM3D street dataset does contain a rough width information under OBJEKTART, but not about the "importance" of the road. Depending on the application, either case could be used.

Possible workflow from a swissTLM3D street shape file with line features to a finished resolved forests shape file, with atmospheric simulations in mind (no emission modeling is intended):

1. Drop all features where the BELAGSART is not paved (e.g. by Select by Expression
2. Drop all features that are underground (STUFE < 0)
3. Create new field RADIUS and assign values based on the OBJEKTART with the field calculator.
4. Buffer the dataset by the RADIUS.
5. If required, assign a field that allows to map a certain amount of traffic to the street type (e.g. classify an actual highway as a highway). The classification mapping is done by the dictionary in the mapdicts-folder.

Crops

relevant namelist flag: docropfields

Manually recorded fields around Yverdon Airfield.

• Required feature type: Polygon
• Required attributes: OBJEKTART
• Optional/Planned attributes: HEIGHT_TOP,CROP,ID

As can be seen in the image in the section Land Cover (Bodenbedeckung), potential land cover datasets do not cover the entire land surface sufficiently to allow for a complete definition of surfaces. This is the case for swisstopo datasets, as well as in Openstreetmap datasets. This is especially often the case where farming is done, as the land cover may quickly change (and corn is significantly different from sugar beet from an aerodynamic and climatic point of view). Originally used for a simulation for Yverdon, where unclassified regions mostly were related to farm land, this enters the whole palmpy ecosystem under the name "crops". As swissTLM3D data also does not include farmland information, a way was found to only include farming information into the static driver. Hence, information provided by this file will appear as vegetation_type = 2 (mixed farming) in the static driver.

As mixed farming is a broad term and, as described above, a corn field has very different parameters for LAI and roughness lenght etc. than a sugar beet field. Therefore, crop classes (as listed in the following table as as guidance) can be individually changed by the namelist parameters veg-/pav-/soil-/wat-/bldg-/albedoparchanges. These parameters also allow individual changes of the base land surface parameters already implemented in PALM.

As mentioned, the following table gives some hints how potential fields could be classified in a GIS. This table can be extended and changed arbitrarily, as 1) all classes will appear as vegetation_type = 2 anyway, and 2) each change to each class needs to be individually specified in the mentioned namelist parameters.

OBJEKTART	Crop Type
1000	Bare Field
1001	Grass Field
1002	Low Crops (Such as beet, Potatoes, )
1003	Wheat
1004	Corn
...	...