diff --git a/index.Rmd b/index.Rmd index f273acd..50bda7f 100644 --- a/index.Rmd +++ b/index.Rmd @@ -28,11 +28,11 @@ code[class^="sourceCode bash"]::before { content: "Bash Source"; } # Vector - Raster ## Learning objectives -* Plot spatial vector and raster data -* Transform vector and raster data -* Write and read spatial vector formats (e.g. KML, GML, GeoJSON, shapefile); -* Apply basic operations on vector data, such as masking, cropping -* Be able to extract raster data using vector data +* Plot spatial vector and raster data; +* Transform vector and raster data; +* Write and read spatial vector formats (e.g. KML, GML, GeoJSON and Shapefile); +* Apply basic operations on vector data, such as masking and cropping; +* Be able to extract raster data using vector data. ## Introduction @@ -46,7 +46,7 @@ The GDAL library is well-documented (http://gdal.org/), but with a catch for R a Thus, functionality that you commonly find in expensive GIS software is also available within R, using free but very powerful software libraries. Here is [handy 'cheatsheet'](https://github.com/rstudio/cheatsheets/blob/master/sf.pdf) for spatial operations with *sf*. The functions of the *sf* package are prefixed by `st_`, short for 'spatial type'. -The possibilities are huge. In this course we can only scratch the surface with some essentials, which hopefully invite you to experiment further and use them in your research. Details can be found in the book *Applied Spatial Data Analysis with R* and several vignettes authored by Roger Bivand, Edzer Pebesma and Virgilio Gomez-Rubio. +The possibilities are huge. In this course we can only scratch the surface with some essentials, which hopefully invites you to experiment further and use them in your research. Details can be found in the book *Applied Spatial Data Analysis with R* and several vignettes authored by Roger Bivand, Edzer Pebesma and Virgilio Gomez-Rubio. This book can be accessed for free through the following [link](https://link.springer.com/book/10.1007/978-1-4614-7618-4)! ### Raster and vector integration and conversion @@ -101,11 +101,11 @@ Raw raster data do not usually conform to any notion of administrative or geogra Step by step we will: -- Download the Landsat 8 data of Wageningen -- Download and prepare administrative boundary data of the Netherlands -- Download Water area data of Wageningen -- Mask the data to match the boundaries of the city -- Mask the data to exclude water bodies +1. Download the Landsat 8 data of Wageningen; +2. Download and prepare administrative boundary data of the Netherlands; +3. Download water area data of Wageningen; +4. Mask the data to match the boundaries of the city; +5. Mask the data to exclude water bodies. ## Prepare the data @@ -140,10 +140,10 @@ We can start by visualizing the data. Since it is a multispectral image, we can wagLandsat[wagLandsat < 0] <- NA # Band names can be changed here -names(wagLandsat) <- c("band1", "band2", "band3", "band4", "band5", "band6", "band7") +names(wagLandsat) <- c("band1", "band2", "band3", "band4", "band5", "band6", "band7") # Select which bands to assign to the red, green, and blue colour channels -plotRGB(wagLandsat, 5, 4, 3) +plotRGB(wagLandsat, 5, 4, 3) ``` We have chosen to visualize the Landsat image as a false color composite, meaning that the chosen bands do not match the RGB channels. Indeed, we have plotted the near-infrared band as red, the red as green, and the green as blue. @@ -207,7 +207,7 @@ In the figure above, the left panel displays the output of `crop`, while the sec We also have a water mask of Wageningen in vector format. Let's download it and also reproject it to the CRS of the Landsat data. ```{block type="alert alert-info"} -Important functions are `st_read` and `st_write`. These are very powerful functions that enable reading and writing simple features or layers from a file or data base. +Important functions are `st_read` and `st_write`. These are very powerful functions that enable reading and writing simple features or layers from a file or database. ``` ```{r} @@ -238,7 +238,7 @@ Also, some of our friends want these exact data too (e.g. the `water` polygon ob One friend of ours is a software engineer and he wants a GeoJSON. Another friend is a GIS-analyst in QGIS and as a backup he wants the file in Geographic Markup Language (GML). These fileformats (GeoJSON and GML, but also KML and Shapefile) are commonly used in spatial analysis. Let’s try to give them the files in those formats! -You can try for yourself and e.g. start by converting them to KML and opening them in Google My Maps (https://mymaps.google.com/). +You can try for yourself and e.g. start by converting them to KML and opening them in [Google My Maps](https://mymaps.google.com/). ```{r, eval=FALSE} diff --git a/index.html b/index.html index 1692754..441e89b 100644 --- a/index.html +++ b/index.html @@ -7,7 +7,7 @@ - + Vector - Raster @@ -24,8 +24,8 @@ } }); -

WUR Geoscripting -WUR logo

+

Vector - Raster

Learning objectives

@@ -988,10 +988,11 @@

Some packages for working with vector and raster data in R

of the sf package are prefixed by st_, short for ‘spatial type’.

The possibilities are huge. In this course we can only scratch the -surface with some essentials, which hopefully invite you to experiment +surface with some essentials, which hopefully invites you to experiment further and use them in your research. Details can be found in the book Applied Spatial Data Analysis with R and several vignettes -authored by Roger Bivand, Edzer Pebesma and Virgilio Gomez-Rubio.

+authored by Roger Bivand, Edzer Pebesma and Virgilio Gomez-Rubio. This +book can be accessed for free through the following link!

Raster and vector integration and conversion

@@ -1199,14 +1200,14 @@

Geometric operations

Visualise spatial vector layer together with a Landsat 8 image from Wageningen

Step by step we will:

- +Netherlands; +
  • Download water area data of Wageningen;
  • +
  • Mask the data to match the boundaries of the city;
  • +
  • Mask the data to exclude water bodies.
  • +

    Prepare the data

    # Install necessary packages and load them
    @@ -1214,9 +1215,9 @@ 

    Prepare the data

    if(!"sf" %in% installed.packages()){install.packages("sf")} library(terra)
    -
    ## terra 1.7.39
    +
    ## terra 1.7.78
    library(sf)
    -
    ## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
    +
    ## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
    # Create data directory if it does not yet exist
     if (!dir.exists("data")) {
       dir.create("data")
    @@ -1237,11 +1238,11 @@ 

    Prepare the data

    wagLandsat[wagLandsat < 0] <- NA # Band names can be changed here -names(wagLandsat) <- c("band1", "band2", "band3", "band4", "band5", "band6", "band7") +names(wagLandsat) <- c("band1", "band2", "band3", "band4", "band5", "band6", "band7") # Select which bands to assign to the red, green, and blue colour channels -plotRGB(wagLandsat, 5, 4, 3)
    -

    +plotRGB(wagLandsat, 5, 4, 3)
    +

    We have chosen to visualize the Landsat image as a false color composite, meaning that the chosen bands do not match the RGB channels. Indeed, we have plotted the near-infrared band as red, the red as green, @@ -1253,8 +1254,7 @@

    Prepare the data

    # Load level 2 dataset nlCitySf <- st_read("data/gadm41_NLD_2.json")
    ## Reading layer `gadm41_NLD_2' from data source 
    -##   `/home/osboxes/Documents/Geoscripting/Tutorials/VectorRaster/data/gadm41_NLD_2.json' 
    -##   using driver `GeoJSON'
    +##   `/home/osboxes/VectorRaster/data/gadm41_NLD_2.json' using driver `GeoJSON'
     ## Simple feature collection with 355 features and 13 fields
     ## Geometry type: MULTIPOLYGON
     ## Dimension:     XY
    @@ -1325,7 +1325,7 @@ 

    Crop, mask and visualise

    # Reset graphical parameters par(opar)
    -

    +

    In the figure above, the left panel displays the output of crop, while the second panel shows the result of masking the Landsat scene using the contour of Wageningen as input.

    @@ -1335,7 +1335,7 @@

    Crop, mask and visualise

    Important functions are st_read and st_write. These are very powerful functions that enable reading and writing simple -features or layers from a file or data base. +features or layers from a file or database.

    # Download and unzip files
    @@ -1344,9 +1344,7 @@ 

    Crop, mask and visualise

    # Load the Water Shapefile directly as an sf object water <- st_read('data/Water.shp')
    -
    ## Reading layer `Water' from data source 
    -##   `/home/osboxes/Documents/Geoscripting/Tutorials/VectorRaster/data/Water.shp' 
    -##   using driver `ESRI Shapefile'
    +
    ## Reading layer `Water' from data source `/home/osboxes/VectorRaster/data/Water.shp' using driver `ESRI Shapefile'
     ## Simple feature collection with 632 features and 32 fields
     ## Geometry type: POLYGON
     ## Dimension:     XY
    @@ -1371,7 +1369,7 @@ 

    Crop, mask and visualise

    # Plot plotRGB(wagLandsatSubW, 5, 4, 3) plot(st_geometry(waterUTM), col = 'lightblue', add = TRUE, border = '#3A9AF0', lwd = 1) # use a site such as https://htmlcolorcodes.com/ to find the hexidecimal code for colours
    -

    +

    Also, some of our friends want these exact data too (e.g. the water polygon object).

    One friend of ours is a software engineer and he wants a GeoJSON. @@ -1380,7 +1378,8 @@

    Crop, mask and visualise

    GML, but also KML and Shapefile) are commonly used in spatial analysis. Let’s try to give them the files in those formats!

    You can try for yourself and e.g. start by converting them to KML and -opening them in Google My Maps (https://mymaps.google.com/).

    +opening them in Google My +Maps.

    # Create output directory if it does not yet exist
     if (!dir.exists("output")) {
       dir.create("output")
    @@ -1437,7 +1436,7 @@ 

    Extract raster values along a transect

    ## max value : 691

    We can start by visualizing the data.

    plot(bel)
    -

    +

    Everything seems correct.

    We want to look at a transect, which we can draw by hand by selecting two points by clicking. The draw('line') function will help @@ -1448,15 +1447,7 @@

    Extract raster values along a transect

    line.

    line <- draw('line')
    ## Loading required package: sp
    -
    ## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
    -## which was just loaded, will retire in October 2023.
    -## Please refer to R-spatial evolution reports for details, especially
    -## https://r-spatial.org/r/2023/05/15/evolution4.html.
    -## It may be desirable to make the sf package available;
    -## package maintainers should consider adding sf to Suggests:.
    -## The sp package is now running under evolution status 2
    -##      (status 2 uses the sf package in place of rgdal)
    -

    +

    Let’s first write this line to a file, using a few different file formats. Note that draw() returns a SpatVector. That means we cannot (directly) use st_write() here like we did @@ -1464,17 +1455,17 @@

    Extract raster values along a transect

    change that using st_as_sf(). After writing to disk, you might try opening these new files in R, or in different software, such as QGIS or ArcGIS.

    -
    # To a GeoJSON
    -st_write(st_as_sf(line), 'output/line.geojson')
    -
    -# To an ESRI Shapefile
    -st_write(st_as_sf(line), 'output/line.shp')
    +
    # To a GeoJSON
    +st_write(st_as_sf(line), 'output/line.geojson')
    +
    +# To an ESRI Shapefile
    +st_write(st_as_sf(line), 'output/line.shp')

    Now the elevation values can simply be extracted using extract().

    -
    alt <- extract(bel, line)
    +
    alt <- extract(bel, line)

    We can plot the result as follows.

    -
    plot(alt$BEL_msk_alt, type = 'l', ylab = "Altitude (m)")
    -

    +
    plot(alt$BEL_msk_alt, type = 'l', ylab = "Altitude (m)")
    +