GSG

Revision as of 14:26, 24 September 2014

A unified framework for environmental monitoring based on a discrete global sampling grid system (GSG)

Background

Environmental monitoring on different spatial scales is a key element in many international and national processes related to ecological and environmental challenges. The increasing availability of large archives of freely available high resolution and Geo-referenced imagery through virtual globes like Google Earth, Microsoft Virtual Earth (bing) or NASA World Wind and others is still relatively underexploited by scientific applications in this context. Even if these information sources are frequently used by scientists, e.g. for communication and visualization of project activities or even as basis for observations in remote sensing imagery, there is hardly any methodological framework that helps to make use of the data sources for sample based estimation. Depending on the available spatial resolution of images, many relevant target variables, like e.g. forest cover, could be interpreted visually and could serve as basis for statistical estimation on larger areas.

Concept of a discrete global sampling grid

We propose a very simple and easy to communicate discrete global grid system (Global Sampling Grid = GSG) that could serve as basic sampling design from global to local scale. The grid construction follows a simple approach: Circles of latitudes are placed with a constant distance on the surface of a sphere with radius of 6378.137 m (major axis of the World Geodetic System, WGS84) outgoing from the equator in North and South direction. On each of these latitudes sampling locations are placed in the same constant distance outgoing from the Greenwich meridian in West and East direction. This leads to a half systematic arrangement of sampling locations, where the distance of points is constant along the latitudes. Here a global GSG with a point distance of 250km is shown in a cylindrical map projection.

The half-systematic design of this grid might become more obvious if draped over a spherical projection of the earth surface, like e.g. used by virtual globes. If exported to a .kml file and shown in Google Earth a GSG with 100km point distance (GSG_100_global) looks like this:

Other examples of GSG subsets can be found on the AWF-Geoserver of the Chair of Forest Inventory and Remote Sensing at the Georg-August-Universität Göttingen. An R package to create grids in any scale and any region of the world has been developed. All grids have the same origin and sampling locations from affine grids of different resolution coincide. This makes the GSG very interesting if data from different sampling studies should be compiled; this is of particular interest for example, for international environmental policy processes and their monitoring instruments, and for internationally acting organizations. A lot of synergies and option for analyzing data are arising if environmental monitoring and sampling studies from other disciplines are combined. One objective of the GSG project is to further develop a R package that allows to construct grids in different spatial resolution and for different regions that are all part of the overall grid system. As preliminary work to the project we developed a basic GSG R code that will be further extended.

GSG application

Given a unified framework to generate sampling grids, the question remains how to make use of this for a sample based estimation of specific target variables? The GSG generates dimensionless sampling locations (points), but how to come to meaningful observations as basis for statistical estimations? A central aspect of the GSG project is to answer this question, depending on different target variables. Contrary to the most common approach of creating model-based wall-to-wall predictions from remote sensing imagery, we are focusing on visual interpretation of relatively small samples. These could be obtained by classifying the dimensionless sampling locations according to a suitable land-cover/land-use classification key that needs to be developed under consideration of the actual image resolution. A simple example would be the estimation of forest cover by the GSG points as observation units (see Estimating forest area). As preliminary work to the project we implemented a small sampling study based on a time series of Landsat lookup images of the district Buritis in Northern Randônia, Brazil. The classification of points from a GSG with 10km resolution into forest or non-forest is basis to estimate forest cover and its changes over time. The respective GSG_10 subset and the imagery used can be found on the AWF-Geoserver (open map and zoom to Brazil). The following figure shows some of the available aquisitions:

An estimate of forest cover change and its standard error over time are given in Figure 1.

estimated forest cover and resulting deforestation rates in Buritis, Northern Randônia, Brazil. Error bars represent the absolute standard error of the estimated cover percentage.

Optimazation of observation units

GSG R code

#############################################################################
#                            GLOBAL SAMPLING GRID  (GSG)                    #
#############################################################################
#by Lutz Fehrmann
 
# Specify your output location and file name for kml export
output.location <- "D:/Projects/GSG/GSG_export/"
file.name <- "test" # e.g. specification of the aoi
 
# Set distance between points in km
d = 100
 
# Optional: Define area of interest (aoi) in a rectangular envelope
# c(min Lat, max Lat, min Lon, max Lon) in decimal degree!
# To compute a global grid use default c(-90, 90, -180, 180)
aoi <- c(-90, 90, -180, 180) 
 
 
# Optional: Compute grid for a specific country: Enter ISO3 country code 
# (http://en.wikipedia.org/wiki/ISO_3166-1_alpha-3#Current_codes)
country <- "CAN"
 
# Calculate grid coordinates ------------
# The grid is systematic (equidistant points along cicles of latitude)
# on a spheroid (WGS84/Pseudo-Mercator, epsg:3857)
# Calculate the angle of an d-arc in degree
deg <- (d/(pi * 6378.137)) * 180 
 
# Create a vector of longitue along equator
longitude <- sort(c(seq(0, -180, -deg), seq(deg, 180, deg)))
 
# Create a vector of latitudes in aoi
latitudes <- longitude[longitude >= aoi[1] & longitude <= aoi[2]]
 
# Calculate a matrix of longitudes for each circle of latitude
longitudes <- outer(longitude, latitudes, function(x,y) x/cos(pi/180 * y))
 
# Compile coordinates 
coordinates<-cbind(as.vector(t(longitudes)), latitudes)
colnames(coordinates)<-c("X","Y")
 
# Subset coordinates in aoi
coordinates<-coordinates[coordinates[,1] >= aoi[3] & coordinates[,1] <= aoi[4],]
 
# To SpatialPointDataFrame
library("rgdal")
GSGsp <- SpatialPointsDataFrame(cbind(coordinates[,1], coordinates[,2]), 
                                data = data.frame(1:nrow(coordinates)))
# Assign projection (longlat +datum=WGS84)
proj4string(GSGsp) <- CRS("+init=epsg:4326") 
 
# Subset gridpoints falling on land (or in a specific country)
library("maptools")
library("rgeos")
 
# load country polygons (from package maptools)
data(wrld_simpl) 
 
# Set projection
wrld_trans <- spTransform(wrld_simpl, CRS("+init=epsg:4326")) 
 
# Subset points on land
GSGland <- GSGsp[wrld_trans,] 
 
#subset per country
GSGcountry <- GSGland[wrld_trans[wrld_trans$ISO3==country,],] 
 
# Plot preview map of aoi ------------
library("maptools")
library("rgeos")
plot(wrld_trans[wrld_trans$ISO3==country,], col = gray(.95), axes = T)
points(GSGcountry, pch = 20, cex = 0.1, col = "red")
 
# Optional: Output of grid coordinates as .kml file for Google Earth or GIS
# Export kml file of the whole global grid with specified point distance -------------
# writeOGR(GSGsp, paste(output.location, "GSG", d, "_global_", file.name, ".kml", sep = ""), layer = paste("GSG_", d, "_kml", sep = ""), driver = "KML")
 
# Export .kml file of the subset of gridpoints on land ------------
# writeOGR(GSGland, paste(output.location, "GSG", d, "_land_", file.name, ".kml", sep = ""), layer = paste("GSG_", d, "_kml", sep = ""), driver = "KML")
 
# Export .kml file of gridpoints for selected country only ----------
writeOGR(GSGcountry, paste(output.location, "GSG", d, "_", country, "_", file.name, ".kml", sep = ""), layer = paste("GSG_", d, "_kml", sep = ""), driver = "KML")