Extracting land use values#
This is a summary of the methods used to extract topographical features from vector and points layers. The source for the vector layer is swissTLMRegio.
QGIS is not the optimal tool for this operation if the intention is to automate the aquisition of the data. However for a limited number of users and relatively few sample locations added each year QGIS is adequate to the task. A summary of the methods for using QGIS is below.
Note: 1’500 meter buffer (circular)
The land-use chapter in the federal report generated alot of interest Land use profile - german. It inspired an academic article and collaboration with Wagenigen Research and University.
Why don’t the percentages add up to 100 ?
This is because some features listed on the map are superimposed on to others. For example a wooded area may contain a public park, this would mean that the park is considered a feature and has a definite surface area but this is in addition to the surface area of the forested area where the park is located. In general we can define the difference in this way:
Land-cover: Are elements such as buildings, forest, orchard or undefined. There is at least one land-cover element for each survey location. If so then it takes 100% of the available dry land.
Land-use: Are elements that may or may not be in the buffer. These are items like schools, hospitals, water-treatment plants. Individually they are only a small portion of the available dry land.
Technical requirements#
In addition to pandas, python and numpy, the following python packages are required:
Geopandas
Shapely
A point and a buffer#
The first step is to create a buffer of the appropriate radius with the point (survey locatio) at the center. The point is constructed by using gps coordinates for the location as arguments to shapely.geometry.Point
. This is then converted to Geopandas.GeodataFrame (gdf) using the coordinate reference system (crs) of the point data. From the gdf, use the buffer method and collect the bounds.
Making a buffer from a point
# method to make the buffer, the point and the bufferbounds
def a_point_buffer(this_location: pd.Series, buffer_radius: int = 1500, coordinate_crs: str = 'epsg:4326', to_crs: str = "epsg:2056"):
# create a point, a buffer and buffer bounds for one location
this_geometry = [Point(this_location.longitude, this_location.latitude)]
this_shape = gpd.GeoDataFrame(this_location, geometry=this_geometry, crs=coordinate_crs)
this_shape = this_shape.to_crs(to_crs)
this_buffer = gpd.GeoDataFrame(geometry=this_shape.buffer(buffer_radius), crs=to_crs)
location_buffer = this_buffer.dissolve()
return this_shape, location_buffer, location_buffer.total_bounds
# select a location
a_location = river_locations[5]
# retrieve the grid coordinates
this_location = locations_and_coordinates.loc[[a_location]].copy()
# make the point, the buffer and the bounds
this_point, abuffer, bufferbounds = a_point_buffer(this_location)
Retrieve location attributes: land-cover#
The bounds and the buffer from the previous operation are used to construct an overlay for the .SHP file that holds the souce data. The name of the location is used to identify the buffer.
def location_geo_attributes(bounds, buffer, shape_file, objcolumn: str = 'OBJVAL', objval: str ='land-cover', buffer_name: str = None):
minx, miny, maxx, maxy = bounds
landcover = gpd.read_file(shape_file, bbox=(minx, miny, maxx, maxy))
landcover = landcover.drop_duplicates()
landcover = gpd.overlay(landcover, buffer, how='intersection', keep_geom_type=True)
landcover.rename(columns={objcolumn:objval}, inplace=True)
landcover['buffer-name'] = buffer_name
return landcover[['buffer-name', objval, 'geometry']]
# indentify the contents of the buffer
buffer_content = location_geo_attributes(bufferbounds, abuffer, 'data/ignorethis/shapes/landcover.shp', buffer_name=a_location)
Undefined areas#
Not all areas of the buffer will have a defined land-cover. To identify those areas, all of the buffer_content
is joined. For this we use the shapely.unary_union
method and provide the buffer_content
as an argument. This will produce a separate gdf that is the negative of the buffer_content
.
# identify the undefined areas
combined_polygon = unary_union(buffer_content_i.geometry)
combined_Pn = gpd.GeoDataFrame({'geometry': [combined_polygon]}, crs=buffer_content_i.crs)
undefined = compute_undefined_sections_of_buffer(abuffer,combined_Pn, buffer_name=a_location)
# combine the known land cover with the undefined
buffer_content_i = pd.concat([buffer_content_i, undefined])
# assign
buffer_content = assign_colors_to_polygons(buffer_content_i, colors=land_use_colors, objval='land-cover')
buffer_content = compute_landuse_percent_of_buffer(buffer_content, abuffer.geometry.area[0])
buffer_content['area_m2'] = buffer_content['area_m2'].astype(int)
To complete the process (for land-cover) the undefined areas are merged with the buffer_content and the surface area for each polygon is calculated.
Public services: land-use#
As noted earlier, some topographical features are overlayed on to other features. Here we are talking about two specific vector layers in the swissTLMRegio collection:
Nutuzungsareal : infrastructure or other usage
Freizeitareal : recreation area
We are interested in these layers because they give us information about the activities that may occur in proximity to the survey location. The same method is used to extract the polygons as was done with land use. However, we are using the shapely.unary_union
method on the entire layer. Which means that we are not differentiating between the subcategories in each layer.
Once the buffer overlay has been applied to both layers the contents are combined and added to the buffer_content data:
# infrastructure
args = dict(objcolumn='OBJEKTART', objval='land-use', buffer_name=a_location)
nz = location_geo_attributes(bufferbounds, abuffer, 'data/ignorethis/shapes/nutuzungsareal.shp', **args)
nz_comb = unary_union(nz.geometry)
nz_comb = gpd.GeoDataFrame({'geometry': [nz_comb]}, crs=buffer_content_i.crs)
# recreational areas
args = dict(objcolumn='OBJEKTART', objval='land-use', buffer_name=a_location)
fz = location_geo_attributes(bufferbounds, abuffer,'data/ignorethis/shapes/freizeit.shp', **args)
fz_comb = unary_union(fz.geometry)
fz_comb = gpd.GeoDataFrame({'geometry': [fz_comb]}, crs=buffer_content_i.crs)
# combine them all into one category called public services
public_service = pd.concat([nz, fz])
ps_comb = unary_union(public_service.geometry)
args = {
'geometry': [ps_comb],
'land-cover':['public-services'],
'buffer-name':[a_location],
'color':'red'
}
ps_comb = gpd.GeoDataFrame(args, crs=buffer_content_i.crs)
ps_comb = compute_landuse_percent_of_buffer(ps_comb, abuffer.geometry.area[0])
# combine the public services and the land-cover
land_use_profile = pd.concat([buffer_content, ps_comb])
land_use_profile['area_m2'] = land_use_profile['area_m2'].astype(int)
land_use_profile['rate'] = land_use_profile['rate'].round(2)
land_use_profile.rename(columns={'rate': '% of buffer'}, inplace=True)
Which results in the land-use profile for the location in question, (not including streets).
Land use and cover for a location.
buffer-name | land-cover | area_m2 | % of buffer |
---|---|---|---|
aare_bern_scheurerk | Wald | 1'387'333 | 0,20 |
aare_bern_scheurerk | Wald | 141'644 | 0,02 |
aare_bern_scheurerk | Siedl | 408'131 | 0,06 |
aare_bern_scheurerk | Siedl | 10 | 0,00 |
aare_bern_scheurerk | Siedl | 75'551 | 0,01 |
aare_bern_scheurerk | Siedl | 671'792 | 0,10 |
aare_bern_scheurerk | Siedl | 2'983'822 | 0,42 |
aare_bern_scheurerk | Siedl | 6'485 | 0,00 |
aare_bern_scheurerk | undefined | 1'382'463 | 0,20 |
aare_bern_scheurerk | public-services | 1'484'991 | 0,21 |
Building a map of land cover#
![_images/bbba01187714460b819dd42491a1bd8fb9a8e916220e4b9639662da6b4548205.png](_images/bbba01187714460b819dd42491a1bd8fb9a8e916220e4b9639662da6b4548205.png)
![_images/a64eb1697d97a2c97f55111174701188cb2264cb2b443b70b9326f404ef577c7.png](_images/a64eb1697d97a2c97f55111174701188cb2264cb2b443b70b9326f404ef577c7.png)
![_images/d79ea3f6f56d33e354b45d14afe3966d8e79b78559816d442d8e4f042e92d6d3.png](_images/d79ea3f6f56d33e354b45d14afe3966d8e79b78559816d442d8e4f042e92d6d3.png)
![_images/caa5de4c89c777eabed78495be6eb2c235343656f67c396dd0e0ba3234b8b69f.png](_images/caa5de4c89c777eabed78495be6eb2c235343656f67c396dd0e0ba3234b8b69f.png)
![_images/5ec1138bd0b8452cb045bc55fe493fefcb44383ddd5134f25d35ef70d56f3e7c.png](_images/5ec1138bd0b8452cb045bc55fe493fefcb44383ddd5134f25d35ef70d56f3e7c.png)
![_images/3ff196b4b5e9be51e3a43e0f0ba8eb3959eeca55706c76d97505c4e93f79809e.png](_images/3ff196b4b5e9be51e3a43e0f0ba8eb3959eeca55706c76d97505c4e93f79809e.png)
![_images/92a28e01dbeb54ecd0bf59c62eb0a5a3b329b6a3c4dee588ec946d95b95c71e9.png](_images/92a28e01dbeb54ecd0bf59c62eb0a5a3b329b6a3c4dee588ec946d95b95c71e9.png)
Using QGIS#
For this method we are using the land-cover layer from swissTLM regio
finished columns = slug, attribute , attribute_type, area, dry, scale
In QGIS:
create a buffer around each survey point
make sure that the survey location and feature_type is in the attributes of the new buffer layer
the survey locations are loaded as points from .csv file
reproject the points layer to the project CRS
use the new buffer layer as an overlay to the land-cover layer
use the overlay intersection tool
select the fields to keep from the buffer (slug and feature type)
select the fields to keep from the land-cover layer
run the function
this creates a temporary layer called intersection
get the surface area of all the land-cover and land-use features in each buffer of the temporary layer
use the field calculator for the attribute table of the layer
in the field calculator, make a new field and enter the formula
\$area
for this example the method is elipsoid bessel 1841 (epsg 7001)
this is set in the properties of the QGIS project
Export the layer as .csv
verify the land-use features per location
drop duplicate values: use location, feature and area to define duplicates
attention! different names for lake and reservoir
change Stausee to See
make a dry land feature
this is the surface area of the buffer that is not covered by water
substract the area of See from the area of the buffer
identify survey locations that have siginifcant water features but are not listed as lakes
Scale the land-use attributes of interest to the available dry-land
Example making dry land columns and scaling the land-use
# separate locations that are lakes
# recall that feature type has a designator for lakes
lakes = lg[(lg.feature_ty == 'l') | lg.slug.isin(snl)].copy()
# from this subset of data separate the surface area covered by water
# set the slug to the index and substract the surface area of the water
# from the surface area of the buffer
lake_only = lakes[lakes.feature == "See"]
lo = lake_only[["slug", "area"]].set_index("slug")
# substract the lake value from the area of the buffer
lo["dry"] = 7065000 - lo.area
lodry = lo["dry"]
# merge the original land-use data from lakes with the
# the dry land results
lgi = lakes.merge(lo["dry"], left_on="slug", right_index=True)
# remove the lake feature from the features columns
lgi = lgi[lgi.feature != "See"].copy()
# scale the landuse feature to the available dry land
lgi["scale"] = (lgi.area/lgi.dry).round(3)
# repeat the process for locations that do not have a lake feature
# these locations are accounted for above
eliminate = [*snl, *lo.index]
# recuperate all other locations
rivers_parcs = lg[~lg.slug.isin(eliminate)].copy()
# define the dry land as the area of the buffer
rivers_parcs["dry"] = 7065000
# scale the features with the dry land value
rivers_parcs["scale"] = rivers_parcs.area/rivers_parcs.dry
# combine the two results
land_cover = pd.concat([rivers_parcs, lgi])
Once the dry land value is calculated for each buffer of the land-cover buffer use the dry-land value to scale the components of the land-use buffer
Extracting street lengths#
The tlmREgio streets or strasse layer is applicable.
Note: the street lenghts are not scaled
finished columns = slug, attribute, attribute_type, length
Attention:
The streets layer can have many copies of the same street in a buffer
The cumulative sum of a street in a layer can exceed the radius of the buffer (windy roads)
Dissolving the polylines of each buffer in the layer is essential
Drop duplicate values on slug, attribute, length
Using the same 1’500 m buffer that was used for land-cover and land-use. Intersect the buffer with the street layer
Dissolve the polylines in each buffer
select the fields from the streets layer to keep (OBJVAL)
select the fields from the buffer layer to kepp (slug, feature_ty)
check the box
keep disjoint features separate
run the function
export to .csv
Land-cover#
Land-cover attributes
array(['undefined', 'Siedl', 'Wald', 'Reben', 'Obstanlage'], dtype=object)
land-cover attribute translations
land_cover_fr = {
'undefined': 'Non défini',
'Siedl': 'Siedl',
'Wald': 'Forêt',
'Reben': 'Vignes',
'Obstanlage': 'Verger'
}
land_cover_en = {
'undefined': 'Undefined',
'Siedl': 'Settlement',
'Wald': 'Forest',
'Reben': 'Vines',
'Obstanlage': 'Orchard'
}
Land-cover results for one location
slug | attribute | attribute_type | area | dry | scale | |
---|---|---|---|---|---|---|
141 | parc-des-pierrettes | undefined | land-cover | 68285 | 3832559 | 0.018 |
623 | parc-des-pierrettes | Siedl | land-cover | 3606885 | 3832559 | 0.941 |
624 | parc-des-pierrettes | Wald | land-cover | 157389 | 3832559 | 0.041 |
Land-use#
Land-use attributes:
array(['Baumschule', 'Friedhof', 'Schul- und Hochschulareal',
'Wald nicht bestockt', 'Abwasserreinigungsareal',
'Historisches Areal', 'Kraftwerkareal', 'Schrebergartenareal',
'Truppenuebungsplatz', 'Unterwerkareal',
'Kehrichtverbrennungsareal', 'Spitalareal',
'Oeffentliches Parkareal', 'Messeareal',
'Massnahmenvollzugsanstaltsareal', 'Kiesabbauareal',
'Steinbruchareal', 'Klosterareal', 'Deponieareal', 'Antennenareal',
'Lehmabbauareal'], dtype=object)
Land-use translations and groups
land_use_fr = {
'Baumschule': 'Pépinière',
'Friedhof': 'Cimetière',
'Schul- und Hochschulareal': 'Zone scolaire et universitaire',
'Wald nicht bestockt': 'Forêt non peuplée',
'Abwasserreinigungsareal': 'Zone de traitement des eaux usées',
'Historisches Areal': 'Zone historique',
'Kraftwerkareal': 'Zone de centrale électrique',
'Schrebergartenareal': 'Zone de jardins familiaux',
'Truppenübungsplatz': 'Terrain d\'entraînement militaire',
'Unterwerkareal': 'Zone de sous-station',
'Kehrichtverbrennungsareal': 'Zone d\'incinération des déchets',
'Spitalareal': 'Zone d\'hôpital',
'Öffentliches Parkareal': 'Zone de parc public',
'Messeareal': 'Zone d\'exposition',
'Massnahmenvollzugsanstaltsareal': 'Zone d\'établissement de traitement',
'Kiesabbauareal': 'Zone d\'extraction de gravier',
'Steinbruchareal': 'Zone de carrière',
'Klosterareal': 'Zone de monastère',
'Deponieareal': 'Zone de décharge',
'Antennenareal': 'Zone d\'antennes',
'Lehmabbauareal': 'Zone d\'extraction d\'argile'
}
land_use_en = {
'Baumschule': 'Nursery',
'Friedhof': 'Cemetery',
'Schul- und Hochschulareal': 'School and University Area',
'Wald nicht bestockt': 'Non-stocked Forest',
'Abwasserreinigungsareal': 'Wastewater Treatment Area',
'Historisches Areal': 'Historical Area',
'Kraftwerkareal': 'Power Plant Area',
'Schrebergartenareal': 'Allotment Garden Area',
'Truppenübungsplatz': 'Military Training Ground',
'Unterwerkareal': 'Substation Area',
'Kehrichtverbrennungsareal': 'Waste Incineration Area',
'Spitalareal': 'Hospital Area',
'Öffentliches Parkareal': 'Public Park Area',
'Messeareal': 'Exhibition Area',
'Massnahmenvollzugsanstaltsareal': 'Correctional Facility Area',
'Kiesabbauareal': 'Gravel Extraction Area',
'Steinbruchareal': 'Quarry Area',
'Klosterareal': 'Monastery Area',
'Deponieareal': 'Landfill Area',
'Antennenareal': 'Antenna Area',
'Lehmabbauareal': 'Clay Extraction Area'
}
# land_uses_grouped:
# outdoor non technical use:
lu_non_tech = ['Friedhof', 'Hitorisches Areal', 'Schrebergartenareal', 'Öffentliches Parkareal', 'Messeareal', 'Klosterareal', 'Wald nicht bestockt']
# technical-extraction-incineration
lu_extraction = ['Kiesabbauareal', 'Steinbruchareal', 'Lehmabbauareal',]
# waste-water-treatment-powere
lu_technical = ['Kehrichtverbrennungsareal', 'Deponieareal', 'Deponieareal', 'Abwasserreinigungsareal','Unterwerkareal', 'Antennenareal']
# services:
lu_serives = ['Massnahmenvollzugsanstaltsareal', 'Schul- und Hochschulareal', 'Spitalareal']
Land-use for one location:
slug | attribute | attribute_type | area | dry | scale | |
---|---|---|---|---|---|---|
685 | parc-des-pierrettes | Abwasserreinigungsareal | land-use | 49488 | 3832559 | 0.012913 |
686 | parc-des-pierrettes | Friedhof | land-use | 11325 | 3832559 | 0.002955 |
687 | parc-des-pierrettes | Oeffentliches Parkareal | land-use | 175591 | 3832559 | 0.045816 |
688 | parc-des-pierrettes | Schrebergartenareal | land-use | 17306 | 3832559 | 0.004516 |
689 | parc-des-pierrettes | Schul- und Hochschulareal | land-use | 1189161 | 3832559 | 0.310279 |
Streets#
Streets attributes:
array(['Autostr', 'Fahrstraes', 'Fussweg', 'HauptStrAB6', 'VerbindStr4',
'NebenStr3', 'VerbindStr6', 'Autobahn', 'HauptStrAB4', 'NebenStr6',
'Autob_Ri'], dtype=object)
Streets: translations and groups
streets_fr = {
'Autostr': 'autoroute',
'NebenStr3': 'route secondaire 3',
'NebenStr6': 'route secondaire 6',
'HauptStrAB6': 'route principale 6',
'HauptStrAB4': 'route principale 4',
'Fahrstraes': 'chemin carrossable',
'Fussweg': 'chemin pédestre',
'Autobahn': 'autoroute',
'Autob_Ri': 'autoroute',
'VerbindStr4': 'route de liason 4',
'VerbindStr6': 'route de liason 6',
}
streets_en = {
'Autostr': 'freeway',
'NebenStr3': 'surface streets 3',
'NebenStr6': 'surface streets 3 6',
'HauptStrAB6': 'inter regional 6',
'HauptStrAB4': 'inter regional 4',
'Fahrstraes': 'bridle path',
'Fussweg': 'pedestrian trail',
'Autobahn': 'freeway',
'Autob_Ri': 'freeway',
'VerbindStr4': 'intra regional 4',
'VerbindStr6': 'intra regional 6',
}
The streets at one survey location:
slug | attribute_type | attribute | length | |
---|---|---|---|---|
638 | parc-des-pierrettes | streets | Autobahn | 1685 |
639 | parc-des-pierrettes | streets | Fahrstraes | 1046 |
640 | parc-des-pierrettes | streets | Fussweg | 3257 |
641 | parc-des-pierrettes | streets | HauptStrAB6 | 2918 |
642 | parc-des-pierrettes | streets | NebenStr3 | 2850 |
643 | parc-des-pierrettes | streets | VerbindStr4 | 67 |
644 | parc-des-pierrettes | streets | VerbindStr6 | 2864 |
Rivers: intersection, size and class#
This requires the 1’500 m buffer created at the begininng and the rivers layer and the lakes layer.
intersection, size and class#
Create a points layer of the intersection of rivers to lakes
Keep only points that are labeled as river (fluss, fluss_u)
Use the attribute table to select the points
Overlay the 1 500 m buffer on the new layer and collect all the intersections within each buffer
export to .csv
Check for intersections that do not have lake name
not all lakes and rivers have names in the layer
find the correct name for the lake and add it to the record
replace river names with nan value with “no name”
distance to intersection#
Note: the intersection points layer needs to be available
Using the
hub lines
tool in QGIS draw lines from survey location to intersect locationLimit the distance to the radius of the buffer
join on the slug field in both layers
eliminate duplicates
calculate the length of the line with
$length
in field calculator of the attribute tableexport to .csv
Identify locations (slugs) that are missing either the river_name or lake
use the previous results (intersection, size and class) to fill in the missing values
check that slugs not in the distance matrix are truly without intersections
slug | lake | river_name | size | class | distance | attribute_type | |
---|---|---|---|---|---|---|---|
72 | parc-des-pierrettes | Le Léman | La Mèbre | 8 | 8 | 333 | river intersection |
Author: hammerdirt-analyst
conda environment: cantonal_report
matplotlib: 3.8.4
pandas : 2.2.2
geopandas : 0.14.4