To install and discover the package, please read the lczexplore_english vignette. To load it simply type:
library(lczexplore)
Usecase: comparing the LCZ available on the WUDAPT european LCZ map and the Geoclimate approach using OpenStreetMap input data
This vignette illustrates the use of the package to compare maps produced with different approaches, it does not aim at drawing conclusion on the relevance of these methods, nor does it aim at describing them in details.
The World Urban Database and Access Portal Tools is a project which aims at creating LCZ with earth observation data. The details of the project and the related papers can be found on the website of the project: https://www.wudapt.org/
The map we are going to use in this vignette is the European LCZ map, as described in the paper available at the following url: https://doi.org/10.1371/journal.pone.0214474
As this method is one of the most used we included the map tif in the embedded data of the package.
The first step is to define a bounding box of the territory of interest. The lczexplore comes with vector data for the Redon Territory, a rural town of Brittany, France.
The importLCZvect is a generic import function for vector data, and when the ouput
argument is set to bBox
, it returns the bounding box of the input layer.
Any bounding box of the bounding box class of the package sf
can be used to crop the raster map. We can use a bounding box created from a vector file. Or create a bounding box from its coordinates. The package osmdata also offers the getbb() function which allow to get the bounding box from a location if it is named in openStreetMap.
redonBbox<-importLCZvect(dirPath=paste0(system.file("extdata", package = "lczexplore"),
"/bdtopo_2_2/Redon"), file="rsu_lcz.geojson",column="LCZ_PRIMARY",output="bBox")
# Example with osmdata function getbb, not run, to avoid installation of the package
# library(osmdata)
# redonBbox<-getbb("Redon", format_out = "sf_polygon")
The second step is to import the raster map and use the previous bounding box to crop the territory of interest. The importLCZraster function performs this in a simple call:
redonWudapt<-importLCZraster(system.file("extdata", package = "lczexplore"),
fileName="redonWudapt.tif", bBox=redonBbox)
#> [1] "/tmp/RtmpOcck5X/Rinstef0a33b0f7ee/lczexplore/extdata/redonWudapt.tif"
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Classes 'sf' and 'data.frame': 2965 obs. of 2 variables:
#> $ EU_LCZ_map: num 12 14 12 12 9 9 12 12 14 14 ...
#> $ geometry :sfc_POLYGON of length 2965; first list element: List of 1
#> ..$ : num [1:5, 1:2] 3415300 3415300 3415231 3415241 3415300 ...
#> ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
#> - attr(*, "sf_column")= chr "geometry"
#> - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
#> ..- attr(*, "names")= chr "EU_LCZ_map"
#> NULL
#> 1 2 3 4 5 6 7 8 9 10 101 102 103 104 105 106 107
The resulting object is an sf file, the pixels of the geotiff have been transformed into geometries. The input file must be a geotiff, whose CRS will be read by the function.
#> Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
#> Coordinate Reference System:
#> User input: ETRS89-extended / LAEA Europe
#> wkt:
#> PROJCRS["ETRS89-extended / LAEA Europe",
#> BASEGEOGCRS["ETRS89",
#> ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
#> MEMBER["European Terrestrial Reference Frame 1989"],
#> MEMBER["European Terrestrial Reference Frame 1990"],
#> MEMBER["European Terrestrial Reference Frame 1991"],
#> MEMBER["European Terrestrial Reference Frame 1992"],
#> MEMBER["European Terrestrial Reference Frame 1993"],
#> MEMBER["European Terrestrial Reference Frame 1994"],
#> MEMBER["European Terrestrial Reference Frame 1996"],
#> MEMBER["European Terrestrial Reference Frame 1997"],
#> MEMBER["European Terrestrial Reference Frame 2000"],
#> MEMBER["European Terrestrial Reference Frame 2005"],
#> MEMBER["European Terrestrial Reference Frame 2014"],
#> ELLIPSOID["GRS 1980",6378137,298.257222101,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[0.1]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4258]],
#> CONVERSION["Europe Equal Area 2001",
#> METHOD["Lambert Azimuthal Equal Area",
#> ID["EPSG",9820]],
#> PARAMETER["Latitude of natural origin",52,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8801]],
#> PARAMETER["Longitude of natural origin",10,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8802]],
#> PARAMETER["False easting",4321000,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8806]],
#> PARAMETER["False northing",3210000,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8807]]],
#> CS[Cartesian,2],
#> AXIS["northing (Y)",north,
#> ORDER[1],
#> LENGTHUNIT["metre",1]],
#> AXIS["easting (X)",east,
#> ORDER[2],
#> LENGTHUNIT["metre",1]],
#> USAGE[
#> SCOPE["Statistical analysis."],
#> AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
#> BBOX[24.6,-35.58,84.73,44.83]],
#> ID["EPSG",3035]]
The LCZ are expressed with several conventions. We chose to use 1 to 10 for the urban and 101 to 107 for the land-use types.
One can then use standard R functions to explore the data.
summary(redonWudapt)
#> EU_LCZ_map geometry geomID
#> 104 :1116 POLYGON :2965 Length:2965
#> 6 : 776 epsg:3035 : 0 Class :character
#> 102 : 701 +proj=laea...: 0 Mode :character
#> 101 : 204
#> 9 : 113
#> 3 : 20
#> (Other): 35
summary(redonWudapt$EU_LCZ_map)
#> 1 2 3 4 5 6 7 8 9 10 101 102 103 104 105 106
#> 0 0 20 0 0 776 0 12 113 0 204 701 0 1116 0 4
#> 107
#> 19
One can see, for instance that no pixel of the imported territory set to one of the following levels: 1, 2, 4, 5, 7, 10, 103 or 105, or that the most present LCZ type is 104, with 1116 pixels.
Once the dataset has been loaded into an R object of class sf
, the function showLCZ
produces the map of it’s repartition. The argumentsf
sets the input dataset, the column of the LCZ types is passed to the argument column
. The argument wf
is just a string that will be used to produce the title of the plot, and it allows the user to specify which workflow produced the input data.
For standard LCZ types, one sets the argument repr="standard"
, and the function will recognize levels from 1 to 10 and from 101 to 107. It is the recommended choice. But standard
representation will also accept levels from 11 to 17 and from A to G.
map1<-showLCZ(sf = redonWudapt,column = 'EU_LCZ_map', repr= "standard")
The following code shows an example of color choice.
These data are embedded in lczexplore
. On can import them using the importLCZvect
function.
This function allows to load geojson or shp files, and one must specify the name of the column containing the LCZ types.
# Path to embedded Data
dirPath<-paste0(system.file("extdata",package="lczexplore"),"/")
# Path to OSM data specifically
dirPathOSM<-paste0(dirPath,"osm/2022/Redon")
# loading Redon Data produced with GeoClimate on OSM input Data
redonOSM<-importLCZvect(dirPath=dirPathOSM,file="rsu_lcz.geojson",column = "LCZ_PRIMARY", drop=TRUE)
summary(redonOSM$LCZ_PRIMARY)
#> 1 2 3 4 5 6 7 8 9 10 101 102 103 104 105 106 107
#> 0 23 2 0 1 33 0 21 60 0 63 2 0 139 138 0 39
The resulting object is an sf object and can therefore also be visualized with the showLCZ function.
map2<-showLCZ(sf = redonOSM, repr="standard")
The purpose of this package is to easily compare local climate zone classifications produced by different algorithms on the same study area.
The compareLCZ
function intersects the geometries of both input sf files. The resulting geometries then either agree or disagree.
Both maps are plotted and a map of where they agree is produced. A confusion matrix is also plotted, revealing how each LCZ type of the first dataset break up into the LCZ types of the second dataset.
If the two datasets are not in the same coordinate reference system (CRS), the function will ask which CRS will be kept.
summary(redonOSM)
#> LCZ_PRIMARY geometry
#> 104 :139 POLYGON :521
#> 105 :138 epsg:32630 : 0
#> 101 : 63 +proj=utm ...: 0
#> 9 : 60
#> 107 : 39
#> 6 : 33
#> (Other): 49
summary(redonWudapt)
#> EU_LCZ_map geometry geomID
#> 104 :1116 POLYGON :2965 Length:2965
#> 6 : 776 epsg:3035 : 0 Class :character
#> 102 : 701 +proj=laea...: 0 Mode :character
#> 101 : 204
#> 9 : 113
#> 3 : 20
#> (Other): 35
compareLCZ(sf1 = redonWudapt, column1 = "EU_LCZ_map", wf1="WUDAPT",
sf2 = redonOSM, column2 = "LCZ_PRIMARY", wf2 = "GeoClimate_OSM",
ref=1 ,exwrite=FALSE,saveG = "comparisonR" )
#> The column EU_LCZ_map of the dataset redonWudapt is the reference against which the LCZ_PRIMARY column of the dataset redonOSM will be compared.
#> Both sf datasets need to live in the same crs projection (srid / epsg),
#> they will be coerced to the specified reference (redonWudapt)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
knitr::include_graphics("comparisonR.png", error=FALSE)
plot of chunk unnamed-chunk-11
The intent of this section is not to discuss the quality of the algorithms used for example, but to illustrate how the package allow to compare grouped categories.
For instance, 51% of the areas of Redon that Wudapt classifies as compact low are seen as compact mid by Geoclimate on OpenStreet map input data. And 34% of what Wudapt classifies as open low are seen as sparsely built by GeoClimate.
But how do Wudapt and GeoClimate agree on broader categories, as built LCZ versus vegetation LCZ ? One could group these LCZs
to explore this, using the groupLCZ
function.
# Regroup levels of the WUDAPT sf file.
redonWudaptGrouped<-groupLCZ(redonWudapt,column="EU_LCZ_map",
urban=c("1","2","3","4","5","6","7","8","9"),
industry="10", vegetation=c("101","102","103","104"),
impervious="105",pervious="106",water="107", outCol = "Regrouped")
# Regroup levels of the GeoClimate sf file
redonOSMgrouped<-groupLCZ(redonOSM, column="LCZ_PRIMARY",
urban=c("1","2","3","4","5","6","7","8","9"),
industry="10", vegetation=c("101","102","103","104"),
impervious="105",pervious="106",water="107", outCol = "Regrouped")
# Pass the files, the columsn and a vector of colors associated to the broader categories
compareLCZ(sf1 = redonWudaptGrouped, column1 = "Regrouped", wf1 = "Wudapt",
sf2 = redonOSMgrouped, column2 = "Regrouped", wf2 = " GeoClimate on OSM data",
exwrite = FALSE, plot=TRUE, repr="alter", levels=c("urban","industry","vegetation","impervious","pervious","water"),colors=c("red","black","green","grey","burlywood","blue"), saveG="comparisonGroupedR")
#> The column Regrouped of the dataset redonWudaptGrouped is the reference against which the Regrouped column of the dataset redonOSMgrouped will be compared.
#> Both sf datasets need to live in the same crs projection (srid / epsg),
#> they will be coerced to the specified reference (redonWudaptGrouped)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
knitr::include_graphics("comparisonGroupedR.png", error=FALSE)
plot of chunk unnamed-chunk-13
One can notice that the two explored algorithms agree rather well about vegetation, as 87% of the surface WUDAPT classifies as vegetation is classified as vegetation by Geoclimate. On urban LCZ, while 68% of the surface WUDAPT classifies as urban is also classified as urban by GeoClimate, 19% is classified as vegetation and 12% as impervious by GeoClimate.
The confusion matrix compare percentage of areas, so if the global area concerned is small, one has to be careful about the interpretation. For instance, the very strong agreement about pervious soils or the relative agreement about water should not be over interpreted as they only represent less than 1% of the areas of the WUDAPT generated map.
For more details on the package, browse other vignettes.
In some case, the source file is a raster with several bands (also called layers).
the function importLCZraster
allows to choose the layer which contains the LCZ values, and when useful
a layer of confidence value. Notice one band was selected by its order in the raster stack, the other by its name, as the function allows both ways.
This allows to also perform a confidence analysis on the agreement between two maps according to the level of confidence one allows to the LCZ value for each pixel. The data used for this example come from Demuzere, M., Kittner, J., Bechtel, B. (2021). LCZ Generator: a web application to create Local Climate Zone maps. Frontiers in Environmental Science 9:637455. https://doi.org/10.3389/fenvs.2021.637455
# Dedicated files were integrated to demonstrate the whole workflow on an area of Sidney
# importation of LCZ map prduced using GeoClimate with OpenStreetMap data
sidneyOSM<-importLCZvect(
dirPath = system.file("extdata/osm/2022/Sidney", package = "lczexplore"), file="sidney_rsu_lcz.geojson",
column="LCZ_PRIMARY",
confid="LCZ_UNIQUENESS_VALUE",
geomID="ID_RSU")
# Creation of the bounding box of this area to crop the raster file produced using the WUDAPT platform
sidneyBbox<-importLCZvect(
system.file("extdata/osm/2022/Sidney", package = "lczexplore"), file="sidney_rsu_lcz.geojson",
,output = "bBox")
# Import of the raster filed, cropped to the bounbding box of the area of interest.
sidneyWudapt<-importLCZraster(system.file("extdata", package = "lczexplore"),
fileName="rasterMD.tif", LCZband=1, LCZcolumn = "LCZraster",
confidenceBand = 3, confidenceColumn = "confidence",
bBox = sidneyBbox )
#> [1] "/tmp/RtmpOcck5X/Rinstef0a33b0f7ee/lczexplore/extdata/rasterMD.tif"
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Classes 'sf' and 'data.frame': 1328 obs. of 3 variables:
#> $ LCZraster : num 5 6 5 3 3 6 6 6 6 3 ...
#> $ confidence: num 64 100 40 92 56 100 100 100 84 56 ...
#> $ geometry :sfc_POLYGON of length 1328; first list element: List of 1
#> ..$ : num [1:5, 1:2] 151 151 151 151 151 ...
#> ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
#> - attr(*, "sf_column")= chr "geometry"
#> - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
#> ..- attr(*, "names")= chr [1:2] "lcz" "classProbability"
#> NULL
#> 1 2 3 4 5 6 7 8 9 10 101 102 103 104 105 106 107
# Make raster confidence value between 0 and 1
sidneyWudapt$confidence<-sidneyWudapt$confidence/100
sidneyComparison<-compareLCZ(
sf1=sidneyOSM, geomID1="ID_RSU", column1="LCZ_PRIMARY", wf1="BD TOPO v2.2", confid1 = "LCZ_UNIQUENESS_VALUE",
sf2=sidneyWudapt, geomID2="geomID", column2="LCZraster", wf2="Open Street Map", confid2="confidence",
ref=1,repr="standard",exwrite=F, location="Part of Sidney")
#> The column LCZ_PRIMARY of the dataset sidneyOSM is the reference against which the LCZraster column of the dataset sidneyWudapt will be compared.
#> Both sf datasets need to live in the same crs projection (srid / epsg),
#> they will be coerced to the specified reference (sidneyOSM)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
confidSidney<-confidSensib(inputDf=sidneyComparison$data, nPoints=5,
wf1="GC on OSM", wf2="WUDAPT",
geomID1="ID_RSU", column1="LCZ_PRIMARY", confid1="LCZ_UNIQUENESS_VALUE",
geomID2="geomID",column2="LCZraster", confid2="confidence",
sep=";", repr="standard",
plot=TRUE, saveG = ".")
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 6.185035e-12 9.896985e-02 2.934190e-01 5.426840e-01 1.000000e+00
#> Warning: Removed 1 rows containing missing values (`geom_point()`).
#> Warning: Removed 1 rows containing missing values (`geom_text()`).
#> Warning: Removed 1 rows containing missing values (`geom_point()`).
#> Warning: Removed 1 rows containing missing values (`geom_text()`).
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.3183997 0.3353561 0.6000000 0.8000000 1.0000000
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.001708454 0.171708619 0.372436250 0.512165544 0.961715868
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.001497977 0.105385470 0.230417890 0.397345488 0.954394752
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.0002142333 0.1103655809 0.3186946386 0.5224514716 0.8800000000
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.0001131941 0.0216191550 0.0453685188 0.1091297317 0.2567207454
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 1.489917e-05 5.954600e-02 1.810297e-01 5.200000e-01 1.000000e+00
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 3.168202e-11 6.386015e-03 2.511515e-02 9.686202e-02 2.104078e-01
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 6.185035e-12 1.046330e-01 1.811952e-01 7.400000e-01 9.808317e-01
#> [1] "confSeq in internFunction "
#> 0% 25% 50% 75% 100%
#> 0.004995776 0.199134045 0.419571234 0.785431982 1.000000000
#> Warning: Removed 9 rows containing missing values (`geom_point()`).
#> Warning: Removed 9 rows containing missing values (`geom_text()`).
#> Warning: Removed 9 rows containing missing values (`geom_point()`).
#> Warning: Removed 9 rows containing missing values (`geom_text()`).
knitr::include_graphics("GeneralUniquenessSensib.png", error=FALSE)
plot of chunk unnamed-chunk-15
One can also explore the evolution of this agreement according to the confidence LCZ type by LCZ type.
knitr::include_graphics("byLCZUniquenessSensib.png", error=FALSE)
plot of chunk unnamed-chunk-16