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:

1. Geopandas

2. 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:

1. Nutuzungsareal : infrastructure or other usage

2. 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

Surface area: 4’145’791 m², % of buffer: 59%

Surface area: 1’528’977 m², % of buffer: 22%

Surface area: 1’484’991 m², % of buffer: 21%

Surface area: 1’382’463 m², % of buffer: 20%

Total length of streets : 17km

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:

1. 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

2. 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

3. 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

4. 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

5. 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

6. 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:

1. The streets layer can have many copies of the same street in a buffer

2. The cumulative sum of a street in a layer can exceed the radius of the buffer (windy roads)

3. Dissolving the polylines of each buffer in the layer is essential

4. Drop duplicate values on slug, attribute, length

5. Using the same 1’500 m buffer that was used for land-cover and land-use. Intersect the buffer with the street layer

6. 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

1. 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

2. Overlay the 1 500 m buffer on the new layer and collect all the intersections within each buffer

   • export to .csv

3. 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

1. Using the hub lines tool in QGIS draw lines from survey location to intersect location

   • Limit 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 table

   • export to .csv

2. 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