palmpy Functions

In this section an overview is given about the functions in the various submodules. Every function implemented in palmpy contains a docstring as well.Therefore, in order to know more about a particular function, it can be examined by help(function) in python with the module imported or by simply looking at the source code.

geodatatools

This submodule contains functions that deal with the processing of geodata. They make heavy use of the gdal library, while more functions could be added using geopandas as well (only one function so far).

• array = cutalti(filein, fileout, xmin, xmax, ymin, ymax, xres, yres)

This function clips raster geodata to specified extents and resolutions and saves on one hand the result in a new file (.tif or .asc), on the other hand it returns the raster as a numpy array.

Example usage:

dhm = cutalti('swissALTI3D2018.tif', 'parentdhm.asc', xmin=720000, xmax=732000, 											ymin=170000, ymax=182000, xres=40, yres=40)

• array = rasterandcuttlm(filein, fileout, xmin, xmax, ymin, ymax, xres, yres, burnatt='OBJEKTART', alltouched=False)

This function rasters a vector shape file, burns in the burnatt attribute and clips and resamples it to the specified exents and resolutions. It saves the raster to a new file (.tif or .asc) and returns the raster as a numpy array. The alltouched flag refers to the rasterization method that takes place. Normally, a raster pixel is only filled with the burn in value if the centroid of the pixel is within the polygon. With alltouched = True all pixels that touch the polygon are filled with the burn value. This makes sense depending on the situation, e.g. if the resolution is chosen as such that streets are not clearly rasterized anymore and a emission simulation is planned. This would significantly alter the results. alltouched = True would ensure that all streets are continously present in the raster, but may be rendered "too wide".

Example usage:

tlmbb = rasterandcuttlm('bb.shp', 'tlm.asc', xmin=730000, xmax=742000, 
                        ymin=190000, ymax=202000, xres=40, yres=40, 											burnatt='OBJEKTART', 'alltouched=False')

• lon, lat = lv95towgs84(E,N)

This function converts LV95 coordinates (or LV03, in which case it adds 2'000'000 and 1'000'000 to the E and N coordinates) to WGS84 coordinates. This is done according to the formula in "Formeln und Konstanten für die Berechnung der Schweizerischen schiefachsigen Zylinderprojektion und der Transformation zwischen Koordinatensystemen" by swisstopo, P.14.

• N, E = epsgtowgs84(E,N)

This function converts any epsg coordinates into WGS84 coordinates and supersedes the above function. This is done using pyproj.

• cutortho(filein, fileout, xmin, xmax, ymin, ymax, xres, yres)

This function clips an orthoimage in geotiff format to the specified extents and resolutions. It does not change a coordinate reference system. This is a basic clipping tool to an external file, without further return objects in python.

• splitroadsshp(path, outpath, majroadlist, minroadlist, attrname = 'OBJEKTART')

This function takes a shape file with street polygons and splits the file into major and minor street files based on the specified attributes and the definitions of major and minor roads (in dictfolder.dialect).

Example usage:

gdt.splitroadsshp(streetsonly, subdir_rasteredshp, mpd.majroads, mpd.minroads)

This step is required, as in a situation, where a minor road is located very close to a highway, it might happen that the less important minor road values are burnt into the raster. Tools to steer this value selection during rasterizations are complex and documentation is sparse.

makestatictools

This submodule contains important definitions and functions for creating a static driver file.

• fillvalues

This dictionary contains the fillvalues for each variable, as PALM requires them to be. Possible fillvalues on the python side are float(-9999.0), int(-9999), np.byte(-127), which are converted to netcdf-compliant NC_FLOAT, NC_INT and NC_FLOAT datatypes when exporting the final static file.

• result = checknestcoordsvalid(dxyp,dxyn,nllx,nlly)

This function checks if the grid spacing of parent and child match. Also it is checked if the lower left position of a nest is aligned with the parent grid (requirement by PALM).

• result = checknxyzvalid(nx,ny,nz)

This function checks if nx, ny and nz values calculated from the specified extents will pass the PALM checks. If not, PALM is not able to start the simulation.

• result = checknxyzisint(nx,ny,nz)

This function checks if the calculated nx, ny and nz values are integer values. Checks for user input errors in the namelist file.

• arr,xll,yll,res = loadascgeodata(filename)

This function loads a raster from a .asc raster file. Usually, when a raster file is processed with a geodatatools function, this function is not necessary. However, when using palmpy without a script, this function lets you import the ascii raster data file previously generated by the geodatatools functions.

Example usage:

topo,xorig,yorig,_ = loadascgeodata('test.asc')

This imports the data from test.asc and saves the altitude information in topo, its origin to xorig and yorig and omits the output of the resolution (by specifying "_"),

• normed, amount = shifttopodown(arr, ischild, shift=None)

This function shifts the topography array downwards, if ischild == 0 then it shifts the topography down so the lowest value reaches 0, if ischild > 0 it shifts the topography down by the supplied shift value. In the latter case, a child only needs to be shifted by the same amount as the parent was shifted downwards.

Example usage:

topo,origin_z = shifttopodown(topo,ischild) 		#when using it for a parent/root
topo,origin_z = shifttopodown(topo,ischild,shift=origin_z) 		#when used for child

• newname = childifyfilename(fileout, ischild)

Modifies the filename to represent its nesting hierarchy (child 1 gets the suffix "_N02" added to its name.)

• vegarr, pavarr, watarr = mapbbclasses2(bbarr, bb2palmvegdict, bb2palmwatdict, fillvalue=[fillvalues['vegetation_type'], fillvalues['water_type']])

This function takes a numpy array (land cover data) and translates the classes to PALM vegetation classes based on bb2palmvegdict and bb2palmwatdict. This function was conceived with swissTLM3D BB data in mind, hence its design to only treat vegetation and water types. However, this approach is in general useful. Pavement arrays are filled with separate routines. Fillvalues are applied if values are not found in the respective dictionaries. This function makes use of the makestatictools.mapdicttoarray() function (added at a later stage, hence the suffix "2" in the function name)

• soilarr = makesoilarray2(vegarr,pavarr, palmveg2soildict, palmpav2soildict, fillvalue = [2,3])

Creates a soi larray from vegetation and pavement array. Everywhere where vegetation_type and pavement_type are nonzero values a soiltype needs to be specified.

This means, that soil types are generically assigned based on the resulting palm vegetation type! If specific soil information is available, read in the shp files with geodatatools, create a new dictionary in mapdicts and perform mapping function. This is not implemented though, as specific soil data has not been come across so far.

The hardcoded fillvalues here specify that, should information be missing in the mapping dictionary, the soil below vegetation will be assigned to class 2 (medium) and for pavement it will be class 3 (medium-fine).

• sfr = makesurffractarray(vegarr,pavarr,watarr)

This function creates a surface_fraction array based on vegetation, pavement and water arrays. Currently in PALM, a pixel can only be either fully vegetation, pavement or water and can therefore only be uniquely defined.

• vegpars,watpars,pavpars,soilpars,bldpars,albpars = createparsarrays(nx,ny)

This function creates vegetation/water/pavement/soil/building/albedo_pars arrays with the correct number of levels and fill values.

index 0: vegetation_pars, 1: water_pars 2: pavement_pars, 3: soil_pars, 4: building_pars, 5:albedo_pars

Example usage to create the soil_pars array:

 _,_,_,soil_pars,_,_ = createparsarrays(nx,ny) 

• parsarr = modifyparsarray(parsarr, npar, newvalue, filterarr, filtervalue)

This function modifies individual parameters on the npar'th level of a _pars-array to a new value depending on a filtervalue in a filterarray.

Example usage:

vegetation_pars = modifyparsarray(vegetation_pars, 10, 12, vegetation_type, 9)

Read: "modify vegetation_pars array at parameter 10 to value 12 where vegetation_type is 9."

• vegpars, albedopars = setalbedovalue(albedopars, vegpars, filterarr, filtervalue, newvalue, npar)

This function combines multiple actions that are required when changing an albedo type. When changing an albedo value of e.g. a vegetation_type, it's albedo_type needs to 1) be set to 0 and 2) the albedo_pars at this pixel needs to be set to the desired albedo_value.

Example usage:

vegpars, albpars = setalbedovalue(albpars, vegpars, vegarr, 9, 0.5, -1)

This changes the vegetation_pars at npar = 10 (albedo type) to 0 where the filtervalue is found in the filterarray, then changes the albedo_pars at these pixels to the desired albedo values. As npar is specified as -1, this is done at all npars in the albedo_pars array.

• newarray = mapdicttoarray(array, mapdict, fillvalues = np.byte(-127))

Maps all entries of an array according to new values given in the mapdict. If key is not found, fillvalue is assigned.

• x, y, nsurface_fraction, nvegetation_pars, nalbedo_pars, npavement_pars, nsoil_pars, nwater_pars, nbuilding_pars = createstaticcoords(xsize, ysize, pixelsize)

Creates static file coordinates that are used to create the xarray DataArrays.

Example usage:

x,_,_,_,_,_,_,_,_ = createstaticcoords(nx[i], ny[i]) 

This creates the x coordinate with correct spacing and values.

• dataarray = createdataarrays(array, dims, coords)

This function turns numpy arrays into xr.DataArrays, that can later be merged to xr.Datasets. To create a DataArray, one needs to specify the numpy data array, the dimensions as a list of strings and the coordinate arrays as a list of variables (in the same order as the dimensions).

Example usage:

vegetation_type = createDataArrays(vegarr,['y','x'],[y,x])

• setneededattributes(dataarray, staticvariable)

Sets the required attributes for each DataArray. Currently, only the minimum set of information is provided and could be extended in the future. staticvariable should be the name of the dataarray (e.g. name zt for topography/zt data).

Example usage:

mst.setneededattributes(zt, 'zt')

• infodict

This dictionary is a basic variant of the global attributes, that will be added to the completely populated Dataset global attributes. The make_static script automatically populates the infodict while running with important data such as origin_x/y/z data.

• setglobalattributes(static, infodict)

This copies data from the infodict to the static Dataset global attributes.

• encodingdict = setupencodingdict(flags)

In order to ensure that each variable is exported to netcdf with the correct datatype, an encoding-dictionary can be supplied. This function creates the correct encoding dict based on which data the user wants to include in the static driver. The flags are set in the namelist.

• outputstaticfile(static, fileout, encodingdict)

This function saves the static xr.Dataset object to a netCDF file in version 3 format with the datatype encodings specified in the encodingdict.

• zlad = createzlad(maxtreeheight, dz)

This function creates the vertical coordinate for the leaf area density array. This coordinate only needs to extend so much to include the hightest trees, which saves a lot of memory. This function is basically the same as createzcoord(), which is described further below. This will function will be deprecated soon.

• showbetadistribution()

This function, when called, will plot various examples of alpha and beta value combinations and their resulting probability density function of the beta distribution. This is intended to be a helper function to select suitable value pairs to create tree crown shapes.

• z = createzcoord(zmax, dz)

This function creates a z coordinate for the static driver. Normally, this is not needed. In order to create a 3D building array though, this is required. Again, this z coordinate only needs to extend so much as to include the highest buildings.

dictfolder.dialect

In the dictfolder, multiple "dialects" can be implemented as separate python files. A dialect in this context is a collection of mapping python dictionaries, that specify which class in the input data format correlates with which PALM class. Which dialect is required can be set in the namelist of the make_static script, which then imports the corresponding file from the dictfolder under the abbreviation mpd. For example, when specifying mapdialect = tlm in the namelist the submodule palmpy.staticcreation.dictfolder.tlm is imported as mpd and each call for mpd will hence load the corresponding mapping dictionary. In this section, the dictionaries that must be specified will be outlined and described.

• bb2palmveg

This dictionary maps the OBJEKTART attribute in the land cover (Bodenbedeckung, BB) shape file to PALM vegetation classes.

If land cover categories have been added manually in the preprocessing stage of a project and new, unique OBJEKTART values have been assigned, they can also be implemented here and mapped to a standard PALM class. It is recommended to assign manually entered classes an unique identifier, such as a number in the thousands (e.g. 1000, 1001 for manually added bare soil or grass areas).

• bb2palmwat

This dictionary maps the OBJEKTART attribute in the land cover (Bodenbedeckung, BB) shape file to PALM water classes. If all this information is already contained in bb2palmveg, this dictionary can be empty.

• pav2palmpav

This dictionary maps the BELAGSART attribute of the assembled pavement shape file to PALM pavement classes.

• str2palmstyp

This dictionary maps the OBJEKTART attribute of the street shape file to PALM street types. The latter is closely linked to the openstreetmap street type classification, but not exactly the same.

• majroads

This list of classes specifies which OBJEKTART values of the street source shape file are regarded as major roads.

• minroads

This list of classes specifies which OBJEKTART values of the street source shape file are regarded as minor roads.

• palmveg2palmsoil

This dictionary governs the assignment of soil type below certain PALM vegetation types.

• palmpav2palmsoil

This dictionary governs the assignment of soil type below certain PALM pavement types.

Adding New Dialects

Currently, there are only a few dialects included in palmpy, it is the tlm and a custom dialect. However, more dialects can be easily implemented by copying and renaming the custom.py file and placing it into the dictfolder location in the palmpy package. To enable the usage of the new dialect, set the mapdialect namelist parameter to the name of the file, without the .py extension. To get the script to recognize your new dialect, look for the following section in the make_static.py script:

NOTE: Currently I would suggest to put your own translation dictionary into the custom.py file. If you have a solid dialect file, let me know and I can implement it into palmpy permanently. Any pull request with changes to dictionaries will be rejected.

if mapdialect == 'tlm':
    import palmpy.staticcreation.dictfolder.tlm as mpd
    print('\nINFO: Imported the tlm mapdict.')
elif mapdialect == 'custom':
    import palmpy.staticcreation.dictfolder.custom as mpd
    print('\nINFO: Imported the custom mapdict.')
elif mapdialect == 'mapdicts':
    import palmpy.staticcreation.mapdicts as mpd
    print('\nINFO: Imported the generic mapdicts mapdict.')
# elif mapdialect == 'NEWDIALECT':
#     import palmpy.staticcreation.dictfolder.NEWDIALECT as mpd
#     print('\nINFO: Imported the NEWDIALECT mapdict.')

To add your newly created dialect to the script, replace the NEWDIALECT entries with the name of your new dialect, uncomment the three lines and you are good to go. Now, the script wil import your created dialect as mpd and as a consequence use your defined dictionaries if you specify it in the namelist.