rosscova rosscova - 2 months ago 9
R Question

Convert map object to array of land/water

I'd like to run some spatial algorithms where it's necessary to "mask out" areas of land. I can make a map showing land:water as black:white with the

maps
package:

lonRange <- c( 100, 150 )
latRange <- c( 10, 50 )
mapObject <- maps::map( database = "world", xlim = lonRange, ylim = latRange, fill = TRUE )


This creates what looks like a fairly simple "list of 4", and displays the data as a black and white map where land is black and water is white. What I'd like is to have that data represented preferably as an array, such that (for example)
land = 0
and
water = 1
. I'm flexible with the exact output, but something like the following would be perfect in my mind. An array where
[,,1]
expresses the binary land/water value,
[,,2]
expresses the latitudes, and
[,,3]
expresses the longitudes:

map <- array( data = NA, dim = c( 10, 10, 3 ) )
map[,,1] <- c( rep( 1, 42 ), rep( 0, 58 ) )
map[,,2] <- c( rep( seq.int( from = latRange[2], to = latRange[1], length.out = 10 ), 10 ) )
map[,,3] <- sort( c( rep( seq.int( from = lonRange[1], to = lonRange[2], length.out = 10 ), 10 ) ) )


So I can call
map[1,1,]
to get the land/sea variable for the given latitude and longitude for that point. In this case 1 (water), at lat_lon 50_100:

> map[1,1,]
[1] 1 50 100


I assume the information I need is contained in the
mapObject
list create above, but I can't figure out how to extract it.

I've been trying to figure out what the values in
mapObject
mean, but can't get my head around them. Elements
[[1]]
and
[[2]]
(
x
and
y
respectively) look promising, but aren't making sense to me enough that I can formulate a plan of how to move forward. The values don't seem to represent latitudes and/or longitudes, and there are
NA
values at the same points in both
x
and
y
.

Can anyone help me decipher the meaning of the contents of
mapObject
? Or suggest a completely different tack to achieve what I'm after here?

The
map
help page says this about the output list:


"the x and y vectors have coordinates of successive polygons"


which should surely help me, but I don't know enough about representations of polygons to know how. If someone can suggest a resource where I could learn what I need here, that would be much appreciated.

Answer

Note: As @hrbrmstr mentioned, there are better methods.

(1) rgeos::gContains() judges whether a point (defined by longitude/latitude) is in SpatialPolygon... (see: Check if point is in spatial object which consists of multiple polygons/holes)(caution: it is a bad idea to handle a lot of ponits in this way). In your case, True means on land. (2) You can make SpatialPolygons from coordinates mapObject has. (Note: mapObject uses NA as a division, so coordinates between NA correspond to one polygon).

library(rgeos); library(maps); library(dplyr)

coords.mat <- matrix(c(mapObject$x, mapObject$y), ncol=2)     # get lon-lat coordinates
div <- c(0, which(is.na(mapObject$x)), length(mapObject$x)+1) # get divisions information

 ## separate coordinates by divisions(NA)
coords.list <- sapply(1:(length(div)-1), function(x) coords.mat[(div[x]+1):(div[x+1]-1),])

 ## change sets of coordinates into SpatialPolygons
map.sp <- coords.list %>% sapply(function(x) Polygon(x)) %>% 
  Polygons(ID = "a") %>% list() %>% SpatialPolygons()

 ## judge
contain.judge <- apply(map, c(1,2), function(x) gContains(map.sp, SpatialPoints(matrix(c(x[3], x[2]), ncol=2))))

 ## chage into binary data and combine
map[,,1] <- as.numeric(contain.judge)

 ## 
plot(coords.mat, pch=".", xlim = lonRange, ylim = latRange)
points(c(map[,,3]), c(map[,,2]), col=c(2,4)[as.factor(c(map[,,1]))], pch = 19)

enter image description here

# bonus: 100x100x3 array version. It takes a few minutes to judge all points.

enter image description here