fieldRS - Tools for remote sensing field work

Why develop fieldRS?

fieldRS¹ 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:

ndvi - RasterStack with 5 time steps of Normalized Difference Vegetation Index (NDVI) images.
fieldData - SpatialPolygonsDataFrame containing ground truth data on crop types.
referenceProfiles - data.frame with NDVI profiles for target crop types.
roads - SpatialLinesDataFrame object with Open Street Map (OSM) street information.

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.

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

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:

diversity - Priority given to the highest class heterogeneity.
richness - Priority given to the highest number of classes.
pixel_frequency - Priority given to the higuest non-NA pixel count.
patch_count - Priority given to the higuest patch count.
road_distance - Priority given to shortest distance.

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
#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

x	y	diversity	richness	patch.count	pixel.frequency	road.distance	ranking
706114.8	4484943	1.398943	1	43	1	1845.245	11
707114.8	4484943	1.437606	1	42	1	1876.150	14
708114.8	4484943	1.418003	1	33	1	1840.659	18
709114.8	4484943	1.347770	1	31	1	1916.704	17
710114.8	4484943	1.442003	1	43	1	2439.187	8

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

x	y	diversity	richness	patch.count	pixel.frequency	road.distance	ranking
706114.8	4484943	1.398943	1	43	1	1845.245	14
707114.8	4484943	1.437606	1	42	1	1876.150	15
708114.8	4484943	1.418003	1	33	1	1840.659	13
709114.8	4484943	1.347770	1	31	1	1916.704	16
710114.8	4484943	1.442003	1	43	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).

#> Warning: Use of `gp$id` is discouraged. Use `id` instead.

#> Warning: Use of `gp$id` is discouraged. Use `id` instead.

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 pixFilter() to apply a erosion filter that removes small clumps likely related to mixed-pixels.

seg.img <- pixFilter(seg.img, 1, "erosion")

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)
#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

The default approach used by extractFields() to draw polygons is the convex hull method, implemented in the chull() function. This is named as the "simple" method in this function. Additionally, the function offers a "complex" approach that applies the concave hull method, implemented in the concaveman() function of the package with the same name. As shown below, the concave hull approach is a bit more precise. However, it is also more time consuming.

fields <- extractFields(seg.img, method="complex")
#> Warning in proj4string(x): CRS object has comment, which is lost in output

#> Warning in proj4string(x): CRS object has comment, which is lost in output

fieldRS - Tools for remote sensing field work

Ruben Remelgado

2020-06-02

Why develop fieldRS?

Example data

Selection of sampling sites

Automated field extraction