fieldRS - Tools for remote sensing field work

Ruben Remelgado


#> Warning: package 'ggplot2' was built under R version 3.5.1
#> Warning: package 'RStoolbox' was built under R version 3.5.1

Why develop fieldRS?

fieldRS1 was designed to ease the collection and management of ground truth data. It provides tools that help select sampling sites, correct and pre-process training data and preview classification results.

Example data

fieldRS contains several raster and vector datasets that can be used to test its tools:

Most data can be accessed with the data() function with the exception of the raster data. This is due to the fact that raster data is stored in the temporary memory when loaded and it can’t be saved as an R object. Below we can see how to load each dataset into R.

data(fieldData) # ground truth data
data(roads) # road shapefile
data(referenceProfiles) # target crop types NDVI profiles
ndvi.ts <- brick(system.file("extdata", "ndvi.tif", package="fieldRS")) # NDVI raster time series

Selection of sampling sites

When aiming for tasks such as land cover classification we try to sample from as many classes as possible in order to build a representative map. However, this is often made difficult by the lack of knowledge of the target area. When selecting sampling sites randomly, this lack of knowledge can lead us to visit homogeneous places that provide us with little information on the composition of the area we wish to classify. To help better select sampling sites we developed rankPlots(). This function uses raster and vector information to prioritize sampling sites based on certain, predefined criteria. As a basis, the function requires a equidistant grid that will define the potential sampling plots. To build this grid we can use the function derivePlots(). This function requires an object from which a spatial extent can be derived (e.g. RasterLayer, SpatialPolygonsDataFrame). Then, based on that extent, the function build a grid of square polygons where its size is predefined by the user. In this example, we will build a grid with 1x1 km polygons. Looking at the output - shown below - we can see that the borders or the image are not considered. This is because those areas result in sampling plots smaller that the target 1x1 km size.

Additionally, rankPlots() requires raster information that provides an indication on the composition of the landscape. For example, the user can choose to use an existence land cover classification. However, when we don’t have existent information at the spatial scale we wish to derive our own land cover map, it can be useful to build an unsupervised classification. To do this, we will use the unsuperCLass() function from the RStoolbox package. We will apply a K-means classifier over our NDVI time series and obtain a cluster image with 5 classes.

k.img <- unsuperClass(ndvi.ts, nSamples=5000, nClasses=5)$map

rankPlots() will use this unsupervised classification to evaluate the potential land cover composition of the landscape. For each potential sampling plot, it will evaluate its fragmentation based on the number of pixel clumps and their size as well as the number of represented classes. As an additional criteria, rankPlots() will allow the inclusion of road information. In practice, the function will evaluate the distance between the closest road and the centre of the sampling plot. In this example, we can use the roads variable provided through fieldRS.

Now that we collected all the reference data, we can use rankPlots() to identify priority sampling sites using the following criteria:

While the function offers a set of predefined sorting criteria, the user can re-arrange them. Let’s show this in action through two test cases. In our first example, we will give priority to the class and patch count. In other words, we will priority the analysis of the unsupervised classification output. In our second example, we will prioritize the distance to the roads. The output will be added to the plot.grid shapefile we previously created.

#> Loading required namespace: igraph
x y class count patch count frequency distance ranking
706114.8 4484943 5 52 1 1845.245 17
707114.8 4484943 5 41 1 1876.150 14
708114.8 4484943 5 30 1 1840.659 18
709114.8 4484943 5 38 1 1916.704 15
710114.8 4484943 5 44 1 2439.187 11
x y class count patch count frequency distance ranking
706114.8 4484943 5 52 1 1845.245 14
707114.8 4484943 5 41 1 1876.150 15
708114.8 4484943 5 30 1 1840.659 13
709114.8 4484943 5 38 1 1916.704 16
710114.8 4484943 5 44 1 2439.187 9

Now let’s plot and compare the output grids coloured by the ranking. Below we can see the output for example 1 (above) and example 2 (below).

Automated field extraction

Once we identify target sampling sites we can fieldRS to select concrete polygons to sample from. For example, if our aim is to map crop types, these polygons might coincide with field parcels. To extract these areas we are going to combine ccLabel()and extractFields(). ccLabel() will allow us to segment a raster image based on the difference between neighbouring pixels. The function will call the mape() function and, for each pixel, will estimate the Mean Absolute Percent Deviation (MAPE) between it and its immediate neighbours using a 3x3 moving window. After iterating through all the pixels in a raster image, the function will apply a user defined threshold to identify breaks between neighbouring field parcels and label each of them using the function clump() from the raster package. In this example, we will use a change threshold of 5%. This means that a pixel will be split from its neighbour if the difference between its value is equal or higher than 5% when compared against the mean estimated between itself and its neighbours. We will perform this analysis over a maximum NDVI composite estimate with ndvi.ts.

ndvi.max <- calc(ndvi.ts, max, na.rm=TRUE) # derive maximum NDVI composite)
seg.img <- ccLabel(ndvi.max, method="spatial", change.threshold=5)$regions # segment NDVI image

Now we can erosionFilter() to remove small clumps that are likely related to mixed-pixels.

seg.img <- erosionFilter(seg.img)

Finally, we will use extractFields() to derive polygons for each segment. This function will draw a polygon based on the extent of each segment. This can be useful when dealing with noisy images as the ones used in this example. However, as the plot below shows, the output might still require some manual editing.

fields <- extractFields(seg.img)