Raster data is like any image. Although it may portray various properties of objects in the real world, these objects don’t exist as separate objects; rather, they are represented using pixels of various values which are assigned a color.
-
Today, it's about constructing a simple spatio-temporal analysis using raster data, R, and git.
+
Today, it’s about constructing a simple spatio-temporal analysis using raster data, R, and git.
Learning objectives
At the end of the lecture, you should be able to:
@@ -895,7 +1038,7 @@
Assumed knowledge from previou
Good scripting habits
Reminder on overall system architecture
-
+
@@ -903,8 +1046,8 @@
Reminder on overall system arch
-
In a previous lecture you got introduced to the overall system architecture of R programming (see figure above) and you saw how to read and write vector data from/to file into your R environment. These vector read/write operations were made possible thanks to the OGR library. The OGR library is interfaced with R thanks to the rgdal package/binding. By analogy, raster data can be read/written thanks to the GDAL library. The figure above provides an overview of the connections between these elements. GDAL stands for Geospatial Data Abstraction Library. You can check the project home page at http://www.gdal.org/ and you will be surprised to see that a lot of the software you have used in the past to read gridded geospatial data use GDAL (i.e.: ArcGIS, QGIS, GRASS, etc). In this lesson, we will use GDAL indirectly via the raster package, which uses rgdal extensively. However, it is possible to call GDAL functionalities directly through the command line from a terminal, which is equivalent to calling a system() command directly from within R. In addition, if you are familiar with R and its string handling utilities, it may facilitate the building of the expressions that have to be passed to GDAL. (Note: This is also doable in Bash scripting, as learned in the previous lesson, and you can even combine the two.)
-Let's start working with rgdal by performing system setup checks.
+
In a previous lecture you got introduced to the overall system architecture of R programming (see figure above) and you saw how to read and write vector data from/to file into your R environment. These vector read/write operations were made possible thanks to the OGR library. The OGR library is interfaced with R thanks to the rgdal package/binding. By analogy, raster data can be read/written thanks to the GDAL library. The figure above provides an overview of the connections between these elements. GDAL stands for Geospatial Data Abstraction Library. You can check the project home page at http://www.gdal.org/ and you will be surprised to see that a lot of the software you have used in the past to read gridded geospatial data use GDAL (i.e.: ArcGIS, QGIS, GRASS, etc). In this lesson, we will use GDAL indirectly via the raster package, which uses rgdal extensively. However, it is possible to call GDAL functionalities directly through the command line from a terminal, which is equivalent to calling a system() command directly from within R. In addition, if you are familiar with R and its string handling utilities, it may facilitate the building of the expressions that have to be passed to GDAL. (Note: This is also doable in Bash scripting, as learned in the previous lesson, and you can even combine the two.)
+Let’s start working with rgdal by performing system setup checks.
The previous function should return the version number of the current version of GDAL installed on your machine. Starting with GDAL 2.0 vector processing becomes incorporated into GDAL. In case the function above returns an error, or if you cannot install rgdal at all, you should verify that all required software and libraries are properly installed. Please refer to the system setup page.
@@ -948,7 +1098,7 @@
Overview of the raster package<
Work on large raster datasets thanks to its built-in block processing functionalities.
Perform fast operations thanks to optimized back-end C code.
Visualize and interact with the data.
-
etc...
+
etc…
@@ -960,13 +1110,13 @@
Overview of the raster package<
-Note: the raster package is now deprecated, as Robert Hijmans has developed a successor to it called terra which is simpler and much faster, as it's rewritten in C++. The downside is that it is not yet supported by the vast majority of other packages, that still assume raster package objects. In the future, this lesson will be rewritten with the use of terra instead of raster.
+Note: the raster package is now deprecated, as Robert Hijmans has developed a successor to it called terra which is simpler and much faster, as it’s rewritten in C++. The downside is that it is not yet supported by the vast majority of other packages, that still assume raster package objects. In the future, this lesson will be rewritten with the use of terra instead of raster.
Explore the raster objects
-
The raster package produces and uses R objects of three different classes. The RasterLayer, the RasterStack and the RasterBrick. A RasterLayer is the equivalent of a single-layer raster, as an R workspace variable. The data themselves, depending on the size of the grid can be loaded in memory or on disk. The same stands for RasterBrick and RasterStack objects, which are the equivalent of multi-layer RasterLayer objects. RasterStack and RasterBrick are very similar, the difference being in the virtual characteristic of the RasterStack. While a RasterBrick has to refer to one multi-layer file or is in itself a multi-layer object with data loaded in memory, a RasterStack may ''virtually'' connect several raster objects written to different files or in memory. Processing will be more efficient for a RasterBrick than for a RasterStack, but RasterStack has the advantage of facilitating pixel based calculations on separate raster layers.
-Let's take a look into the structure of these objects.
+
The raster package produces and uses R objects of three different classes. The RasterLayer, the RasterStack and the RasterBrick. A RasterLayer is the equivalent of a single-layer raster, as an R workspace variable. The data themselves, depending on the size of the grid can be loaded in memory or on disk. The same stands for RasterBrick and RasterStack objects, which are the equivalent of multi-layer RasterLayer objects. RasterStack and RasterBrick are very similar, the difference being in the virtual characteristic of the RasterStack. While a RasterBrick has to refer to one multi-layer file or is in itself a multi-layer object with data loaded in memory, a RasterStack may ‘’virtually’’ connect several raster objects written to different files or in memory. Processing will be more efficient for a RasterBrick than for a RasterStack, but RasterStack has the advantage of facilitating pixel based calculations on separate raster layers.
+Let’s take a look into the structure of these objects.
From the metadata displayed above, we can see that the RasterLayer object contains all the properties that geo-data should have; that is to say a projection, an extent and a pixel resolution.
The RasterBrick metadata displayed above are mostly similar to what we saw earlier for the RasterLayer object, with the exception that these are multi-layer objects.
@@ -1046,7 +1196,7 @@
Reading and writing from/to file
## Unpack the archive
unzip('gewata.zip')
-
Gewata is the name of the data set added, it is a multi-layer GeoTIFF object, its file name is LE71700552001036SGS00_SR_Gewata_INT1U.tif, informing us that this is a subset from a scene acquired by the Landsat 7 sensor. Let's not worry about the region that the data covers for now, we will find a nice way to discover that later on in the tutorial. See the example below.
+
Gewata is the name of the data set added, it is a multi-layer GeoTIFF object, its file name is LE71700552001036SGS00_SR_Gewata_INT1U.tif, informing us that this is a subset from a scene acquired by the Landsat 7 sensor. Let’s not worry about the region that the data covers for now, we will find a nice way to discover that later on in the tutorial. See the example below.
Now that we have downloaded and unpacked the GeoTIFF file, it should be present in our working directory. We can investigate the content of the working directory (or any directory) using the list.files() function.
-
+
+
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
+## prefer_proj): Discarded datum unknown in Proj4 definition
+
-Let's take a look at the structure of this object.
+Let’s take a look at the structure of this object.
+
+
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
+## prefer_proj): Discarded datum unknown in Proj4 definition
+
Note that in addition to supporting most commonly used geodata formats, the raster package has its own format. Saving a file using the .grd extension ('filename.grd') will automatically save the object to the raster package format. This format has some advantages when performing geo processing in R (one advantage for instance is that it conserves original filenames as layer names in multilayer objects), however, it also has disadvantages, since those files are not compressed and thus very large.
-
Geo processing, in memory vs. on disk
-
When looking at the documentation of most functions of the raster package, you will notice that the list of arguments is almost always ended by .... These 'three dots' are called an ellipsis; it means that extra arguments can be passed to the function. Often these arguments are those that can be passed to the writeRaster() function; meaning that most geo-processing functions are able to write their output directly to file, on disk. This reduces the number of steps and is always a good consideration when working with big raster objects that tend to overload the memory if not written directly to file.
+
Note that in addition to supporting most commonly used geodata formats, the raster package has its own format. Saving a file using the .grd extension (‘filename.grd’) will automatically save the object to the raster package format. This format has some advantages when performing geo processing in R (one advantage for instance is that it conserves original filenames as layer names in multilayer objects), however, it also has disadvantages, since those files are not compressed and thus very large.
+
Geo processing, in memory vs. on disk
+
When looking at the documentation of most functions of the raster package, you will notice that the list of arguments is almost always ended by .... These ‘three dots’ are called an ellipsis; it means that extra arguments can be passed to the function. Often these arguments are those that can be passed to the writeRaster() function; meaning that most geo-processing functions are able to write their output directly to file, on disk. This reduces the number of steps and is always a good consideration when working with big raster objects that tend to overload the memory if not written directly to file.
Data type is (still) important
When writing files to disk using writeRaster() or the filename = argument in most raster processing functions, you should set an appropriate data type. Use the datatype = argument, it will save some precious disk space, and increase read and write speed.
See details in ?dataType.
@@ -1156,14 +1317,14 @@
Creating layer stacks
## Retrieve the content of the tura sub-directory
list <- list.files(path='tura/', full.names=TRUE)
-
The object list contains the file names of all the single layers we have to stack. Let's open the first one to visualize it.
+
The object list contains the file names of all the single layers we have to stack. Let’s open the first one to visualize it.
plot(raster(list[1]))
-
We see an NDVI layer, with the clouds masked out. Now let's create the RasterStack, the function for doing that is called stack(). Looking at the help page of the function , you can see that it can accept a list of file names as argument, which is what the object list represents. So we can very simply create the layer stack by running the function.
+
We see an NDVI layer, with the clouds masked out. Now let’s create the RasterStack, the function for doing that is called stack(). Looking at the help page of the function , you can see that it can accept a list of file names as argument, which is what the object list represents. So we can very simply create the layer stack by running the function.
-
Now that we have our 166 layers RasterStack in memory, let's write it to disk using the writeRaster() function. Note that we decide here to save it as .grd file (the native format of the raster package); the reason for that is that this file format conserves original file names (in which information on dates is written) in the individual band names. The data range is comprised between -10000 and +10000, therefore such a file can be stored as signed 2 byte integer (INT2S).
+
Now that we have our 166 layers RasterStack in memory, let’s write it to disk using the writeRaster() function. Note that we decide here to save it as .grd file (the native format of the raster package); the reason for that is that this file format conserves original file names (in which information on dates is written) in the individual band names. The data range is comprised between -10000 and +10000, therefore such a file can be stored as signed 2 byte integer (INT2S).
According to the algorithm description, water is coded as 1, cloud as 4 and cloud shadow as 2.
Does the cloud mask fit with the visual interpretation of the RGB image we plotted before?
-
We can also plot the two on top of each other, but before that we need to assign no values (NA) to the 'clear land pixels' so that they appear transparent on the overlay plot.
+
We can also plot the two on top of each other, but before that we need to assign no values (NA) to the ‘clear land pixels’ so that they appear transparent on the overlay plot.
@@ -1359,7 +1527,7 @@
About the area
## Remove fmask layer from the Landsat stack
tahiti6 <- dropLayer(tahiti, 7)
-
We will first do the masking using simple vector arithmetic, as if tahiti6 and fmask were simple vectors. We want to keep any value with a 'clean land pixel' flag in the cloud mask; or rather, since we are assigning NAs, we want to discard any value of the stack which has a corresponding cloud mask pixel different from 0. This can be done in one line of code.
+
We will first do the masking using simple vector arithmetic, as if tahiti6 and fmask were simple vectors. We want to keep any value with a ‘clean land pixel’ flag in the cloud mask; or rather, since we are assigning NAs, we want to discard any value of the stack which has a corresponding cloud mask pixel different from 0. This can be done in one line of code.
+
+
+# [WUR Geoscripting](https://geoscripting-wur.github.io/) 
