A common task in spatial data analysis is extracting SpatialPoints inside a set of polygons or buffer zones. Analysts can use standard GIS or map tools to extract a set of points within an area of interest using manual “point-and-click” routines. This method is easy, but will probably prove impractical, especially in cases involving big data. The alternative is to train a machine to automatically extract the points in a polygon or buffer zone. This post achieves that task and presents a case-study with R code.
Case Study – Intro
Satellite sensors provide long-term solar radiation and weather data at fixed points across the Earth’s surface. In order to use the satellite data, it is necessary to identify and extract a specific set of SpatialPoints. The task can be especially challenging when:
- There are over 250 million spatial points to choose from, and
- Irregular buffer zones or polygons define which points are relevant.
Fortunately, spatial analysis methods in R can simplify the task. Consider the following case study:
- A solar project developer has found 6 sites in Qatar that are suitable for future use.
- Buffer zones are defined around each site with a 5km and 10km radius.
- In some cases the buffer zones are stand-alone circles. In other cases the buffer zones overlap, creating irregular shaped areas of interest
- The developer must identify SpatialPoints within the buffer zones that represent satellite data collection points.
- Each point contains 20 year time series data on solar radiation and other weather variables essential for estimating solar power production
- The developer must identify all points that fall within the buffer zones. Next, he must extract those points and their meta-data IDs to download the satellite data.
Case Study – Data Definitions and Loading
The first step is to define and download all relevant map data. Map components will include:
- The spatial points for the 6 potential project sites, which are defined manually;
- The lines and polygons that form the base map for the the State of Qatar1
- The spatial points where satellite data collection takes place2
The following R code defines the coordinates for the 6 project sites, loads the base map data and creates 5km and 10km boundaries around the project sites:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 |
# 0.Initialize Code # Load librariers and DLLs library(sp) library(rgdal) library(maptools # paths – update and point to input map vector files data.path <- "/home/bxhorn/Documents/DohaR/data/" # map projection and graph devices map.proj <- CRS("+proj=longlat +ellps=WGS84") op <<- par() ##----------------------------------------------------------------------------------------## # 1.Define spatial points for 6 potential development sites x <- c(51.25, 51.40, 51.25, 51.05, 51.07, 51.05) y <- c(26.00, 25.95, 25.75, 25.30, 25.31, 25.20) coords <- cbind(x, y) site.names <- data.frame(site = LETTERS[1:6]) # go spatial spdf.sites <- SpatialPointsDataFrame(coords, site.names) ##----------------------------------------------------------------------------------------## # 2.Load Qatar map data qatar.bbox <- matrix(c(50.7, 51.65, 24.4, 26.4), ncol = 2, byrow = TRUE, dimnames = list(c("x", "y"), c("min", "max"))) # GSHHS Shores Region Map (SpatialPolygons) qatar.shores <- Rgshhs(paste0(data.path, "gshhs_f.b"), xlim = c(50.65, 51.7), ylim = c(24.4, 26.4)) q.shoreGSHHS <- qatar.shores$SP # NOAA Borders (SpatialLinesDataFrame) WDBIIinfo <- ogrInfo(data.path, "WDBII_border_f_L1") WDBII_b_L1 <- readOGR(data.path, "WDBII_border_f_L1") q.borderNOAA <- pruneMap(WDBII_b_L1, xlim = c(50.65, 51.7), ylim = c(24.4, 26.4)) # GADM Amin Boundries GADMinfo <- ogrInfo(data.path, "QAT_adm1") q.adminGADM <- readOGR(data.path, "QAT_adm1") # Convert from SpatialPolygonDataFrame to SpatialLineDataFrame q.adminGADM <- as(q.adminGADM, "SpatialLinesDataFrame") ##----------------------------------------------------------------------------------------## # 3.Load LSA-SAF pixels point and geographical coordinates lat.mena <- readGDAL(paste0(data.path, "HDF5_LSASAF_MSG_LAT_NAfr_4bytesPrecision")) lon.mena <- readGDAL(paste0(data.path, "HDF5_LSASAF_MSG_LON_NAfr_4bytesPrecision")) # clean-up missing data is.na(lat.mena@data$band1) <- lat.mena@data$band1 == 900000 is.na(lon.mena@data$band1) <- lon.mena@data$band1 == 900000 # define spatial points for satellite pixels over Qatar land mass lon.sp <- lon.mena@data$band1/10000 lat.sp <- lat.mena@data$band1/10000 mat.sp <- cbind(lon.sp, lat.sp) row.names(mat.sp) <- 1:nrow(mat.sp) sp.pixels <- SpatialPoints(na.omit(mat.sp), bbox = qatar.bbox, proj4string = map.proj) ##----------------------------------------------------------------------------------------## 4. Create 10km and 5km Buffer Zones (width in decimal degrees) buff10k <- gBuffer(spdf.sites, width = 0.10) buff5k <- gBuffer(spdf.sites, width = 0.05) |
Case Study – Base Map Creation
The loaded data can next be applied to create a base map of Qatar with the 6 project sites, the satellite pixel points and the buffer zones around each project:
Base map of Qatar with SpatialPoints and buffer zones
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
# 4.Create State of Qatar Map with Satellite Pixels, Proposed Project Sites and Buffer Zones M0title <- "State of Qatar" M1title <- "Project Sites vs. Meteosat-10 Pixel Locations" S0title <- "Pixel Map: MSG Data Servers, Satellite Application Facility, EUMETSAT" S1title <- "Land Map: GSHHS and WDBII Databases, National Geophysical Data Center, NOAA" par(op) par(mar = c(2,2,2,2) + 0.01, pin = c(3.73, 7.5)) plot(q.shoreGSHHS, xaxs = "i", yaxs = "i", axes = FALSE, bg = "lightcyan", col = "papayawhip", border = "lightskyblue3") plot(q.borderNOAA, lwd = 1.3, col = "grey45", add = TRUE) plot(q.adminGADM, lwd = 0.5, col = "grey60", add = TRUE) plot(sp.pixels, col = "plum", cex = 0.45, add = TRUE) plot(spdf.sites, col = "red", pch = 16, cex = 0.95, add = TRUE) plot(buff5k, lty = 2, border = "red", Add = TRUE) plot(buff10k, lty = 2, border = "red", Add = TRUE) axis(1, at = c(50.8 + 0:5 * 0.2), pos = 24.393, las = 1, cex.axis = 0.5, col = "grey20", col.axis = "grey20", labels = parse(text = sp:::degreeLabelsEW(c(50.8 + 0:5 * 0.2)))) axis(2, at = c(24.5 + 0:10 * 0.2), pos = 50.65, las = 2, cex.axis = 0.5, col = "grey20", col.axis = "grey20", labels = parse(text = sp:::degreeLabelsNS(c(24.5 + 0:10 * 0.2)))) axis(3, at = c(50.8 + 0:5 * 0.2), pos = 26.2927, las = 1, cex.axis = 0.5, col = "grey20", col.axis = "grey20", labels = parse(text = sp:::degreeLabelsEW(c(50.8 + 0:5 * 0.2)))) axis(4, at = c(24.5 + 0:10 * 0.2), pos = 51.695, las = 2, cex.axis = 0.5, col = "grey20", col.axis = "grey20", labels = parse(text = sp:::degreeLabelsNS(c(24.5 + 0:10 * 0.2)))) title(main = M0title, font.main = 2, line = 3.25) title(main = M1title, cex.main = 0.85, font.main = 2, line = 2.25) title(sub = S0title, cex.sub = 0.80, font.sub = 3, line = 2.0) title(sub = S1title, cex.sub = 0.80, font.sub = 3, line = 3.00) box(col = "grey20", lwd = 1.3) text(x = 51.1, y = 24.9, label = "Jarayan Al Batnah", cex = .6) text(x = 51.46, y = 25.0, label = "Al Wakrah", cex = .6) text(x = 50.95, y = 25.4, label = "Al Jumayliyah", cex = .6) text(x = 51.15, y = 25.785, label = "Al Ghuwayriyha", cex = .6) text(x = 51.431, y = 25.235, label = "Ar Rayyan", cex = .6) text(x = 51.58, y = 25.32, label = "Ad Dawah", cex = .6) text(x = 51.405, y = 25.46, label = "Umm Salal", cex = .6) text(x = 51.39, y = 25.7, label = "Al Khawr", cex = .6) text(x = 51.21, y = 25.95, label = "Ash Shamal", cex = .6) ##----------------------------------------------------------------------------------------## |
Case Study – Extracting SpatialPoints
And finally, the satellite pixel points within the project buffer zones can be easily extracted as follow:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
# 5. Extract Satellite Pixels Inside the 10km and 5km Buffer Zones # first ensure the required data objects have the same CRS proj4string(sp.pixels) <- WGS84 proj4string(buff10k) <- WGS84 proj4string(buff5k) <- WGS84 # Next two lines confirm the beauty of sp methods # A few dozen satellite pixel points are easily extracted from hundreds of millions!!! ftplist.10km <- sp.pixels[buff10k] ftplist.5k <- sp.pixels[buff5k] ##----------------------------------------------------------------------------------------## # 6. Confirm pixel ID# and coordinates for FTP download ftplist.5k@coords lon.sp lat.sp 674139 51.27 26.03 676350 51.24 26.00 678562 51.28 25.97 678564 51.41 25.97 680775 51.38 25.94 680776 51.44 25.94 682987 51.42 25.91 691830 51.24 25.76 694041 51.22 25.73 694042 51.28 25.73 718364 51.05 25.35 720575 51.03 25.32 720576 51.09 25.32 722787 51.07 25.29 727209 51.01 25.22 727210 51.08 25.22 729421 51.05 25.19 |
The method for extracting is simple and straightforward. Index extraction is used on the satellite SpatialPoints with the [] operator, and within brackets is the name of spatial object for the 5km or 10km buffer zones. Nothing else is required!
- Polygons that define the Qatar shoreline can be downloaded form NOAA at http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html; polygons that define country boundaries can be downloaded from the CIA World Database II at https://www.evl.uic.edu/pape/data/WDB/; and polygons defining the administrative boundaries within the State of Qatar can be downloaded from http://www.gadm.org ↩
- Data that defines the satellite sensor points for the Meteosat-1o satellite can be downloaded from https://landsaf.meteo.pt/auxiliarDataFiles.jsp ↩