Show code cell source Hide code cell source

%load_ext watermark
import pandas as pd
import numpy as np

Extracting land use values

This is a summary of the methods to extract topographical features from vector and points layers using QGIS. Atribute names are translated to the target languages and Null values are either corrected or eliminated.

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 this method is adequate to the task.

Note: 1’500 meter buffer (circular)

The land-use chapter in the federal report generated alot of interest. It inspired an academic article and collaboration with Wagenigen Research and University. The principal was applied using an empirical Bayes method with the Solid-Waste-Team at the EPFL.

1. the area of the buffer is \(\pi * 1500²\) or 7065000 m²

Land-use, land-cover and streets

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.

Extracting land-cover and land-use:

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

# locations with significant still water feature but
# not listed as a lake
see_not_lake = l_u[(l_u.feature == "See")&(l_u.feature_ty != "l")]
snl = see_not_lake.slug.unique()

# lakes and locations in snl
# 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

Show code cell source Hide code cell source

def collect_vitals(data):
    total = data.quantity.sum()
    median = data.pcs_m.median()
    samples = data.loc_date.nunique()
    ncodes = data.code.nunique()
    nlocations = data.slug.nunique()
    nbodies = data.feature_name.nunique()
    ncities = data.city.nunique()
    min_date = data["date"].min()
    max_date = data["date"].max()
    
    return total, median, samples, ncodes, nlocations, nbodies, ncities, min_date, max_date

def find_missing(more_than, less_than):
    return np.setdiff1d(more_than, less_than)
def find_missing_loc_dates(done, dtwo):
    locs_one = done.loc_date.unique()
    locs_two = dtwo.loc_date.unique()
    return find_missing(locs_one, locs_two)

def aggregate_gcaps_gfoams_gfrags(data, codes,columns=["Gfoams", "Gfrags", "Gcaps"]):
    for col in columns:
        change = codes.loc[codes.parent_code == col].index
        data.loc[data.code.isin(change), "code"] = col
        
    return data

def make_a_summary(vitals, add_summary_name=False):

    a_summary = f"""
    Number of objects: {vitals[0]}
    
    Median pieces/meter: {vitals[1]}
    
    Number of samples: {vitals[2]}
    
    Number of unique codes: {vitals[3]}
    
    Number of sample locations: {vitals[4]}
    
    Number of features: {vitals[5]}
    
    Number of cities: {vitals[6]}
    
    Start date: {vitals[7]}
    
    End date: {vitals[8]}
    
    """

    if add_summary_name:
        a_summary = f"""
        Summary name = {add_summary_name}

        {a_summary}
        """
        
    return a_summary
def combine_survey_files(list_of_files):

    files = []
    for afile in list_of_files:
        files.append(pd.read_csv(afile))
    return pd.concat(files)

def indexed_feature_data(file, index: str = "code"):
    df = pd.read_csv(file)
    df.set_index(index, drop=True, inplace=True)
    return df
code_cols = ['material', 'description', 'source', 'parent_code', 'single_use', 'groupname']

group_by_columns = [
    'loc_date', 
    'date', 
    'feature_name', 
    'slug',     
    'parent_boundary',
    'length',
    'groupname',
    'city',
    'code', 
]
agg_this = {
    "quantity":"sum",
    "pcs_m": "sum"
}

survey_data = [
    "data/end_process/after_may_2021.csv",
    "data/end_process/iqaasl.csv",
    "data/end_process/mcbp.csv",
    "data/end_process/slr.csv",
]

code_data =  "data/end_process/codes.csv"
beach_data = "data/end_process/beaches.csv"
land_cover_data = "data/end_process/land_cover.csv"
land_use_data = "data/end_process/land_use.csv"
street_data = "data/end_process/streets.csv"
intersection_attributes = "data/end_process/river_intersect_lakes.csv"
surveys = combine_survey_files(survey_data)
codes = indexed_feature_data(code_data, index="code")
beaches = indexed_feature_data(beach_data, index="slug")
land_cover = pd.read_csv(land_cover_data)
land_use = pd.read_csv(land_use_data)
streets = pd.read_csv(street_data)
river_intersect_lakes = pd.read_csv(intersection_attributes)

	loc_date	date	feature_name	slug	code	pcs_m	quantity	parent_boundary	length	groupname	city	feature_type
0	('tiger-duck-beach', '2021-10-07')	2021-10-07	lac-leman	tiger-duck-beach	G1	0.00	0	rhone	30	food and drink	Saint-Sulpice (VD)	l
1	('tiger-duck-beach', '2021-10-07')	2021-10-07	lac-leman	tiger-duck-beach	G10	0.07	2	rhone	30	food and drink	Saint-Sulpice (VD)	l
2	('tiger-duck-beach', '2021-10-07')	2021-10-07	lac-leman	tiger-duck-beach	G100	0.30	9	rhone	30	waste water	Saint-Sulpice (VD)	l
3	('tiger-duck-beach', '2021-10-07')	2021-10-07	lac-leman	tiger-duck-beach	G101	0.00	0	rhone	30	personal items	Saint-Sulpice (VD)	l
4	('tiger-duck-beach', '2021-10-07')	2021-10-07	lac-leman	tiger-duck-beach	G102	0.00	0	rhone	30	personal items	Saint-Sulpice (VD)	l

Land-cover

Null values:

slug              False
attribute         False
attribute_type    False
area              False
dry               False
scale             False
dtype: bool

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

Null values:

slug              False
attribute         False
attribute_type    False
area              False
dry               False
scale             False
dtype: bool

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

Null values

streets.isna().any()

slug              False
attribute_type    False
attribute         False
length            False
dtype: bool

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',   
}

# surface streets

str_surface = ['NebenStr3', 'NebenStr6']
str_ped_br = ['Fahrstraes', 'Fussweg']
str_main = [ 'HauptStrAB6', 'VerbindStr4','VerbindStr6', 'HauptStrAB4', 'NebenStr6']
str_auto = ['Autobahn', 'Autostr', 'Autob_Ri']

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

slug              False
lake              False
river_name        False
size              False
class             False
distance          False
attribute_type    False
dtype: bool

import reportclass as r_class
import setvariables as conf_
# updating language maps
# Collecting required data to establish a report
# This includes the language maps for all the common
# abbreviations and columns or index labels.
c_l = r_class.language_maps()
# luen = conf_.land_use_fr
# land_cover_fr = {
#     'undefined': 'Non défini',
#     'Siedl': 'Siedl',
#     'Wald': 'Forêt',
#     'Reben': 'Vignes',
#     'Obstanlage': 'Verger'
# }

# luen.update(land_cover_fr)
# fren = pd.read_csv('data/end_process/fr_labels.csv').set_index('en')
# ty = pd.DataFrame(luen.values(), index=luen.keys())
# ty['fr'] = ty[0]
# ty.index.name='en'

# ty.drop(0, axis=1,inplace=True)

c_l['fr'].loc['Siedl']

fr    Siedl
Name: Siedl, dtype: object

# updated = pd.concat([fren, ty], axis=0)
# updated.to_csv('data/end_process/fr_labels.csv', index=True)

%watermark -a hammerdirt-analyst -co --iversions

Author: hammerdirt-analyst

conda environment: cantonal_report

numpy : 1.25.2
pandas: 2.0.3