rosscova rosscova - 1 year ago 78
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


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
water = 1
. I'm flexible with the exact output, but something like the following would be perfect in my mind. An array where
expresses the binary land/water value,
expresses the latitudes, and
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( from = latRange[2], to = latRange[1], length.out = 10 ), 10 ) )
map[,,3] <- sort( c( rep( from = lonRange[1], to = lonRange[2], length.out = 10 ), 10 ) ) )

So I can call
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
list create above, but I can't figure out how to extract it.

I've been trying to figure out what the values in
mean, but can't get my head around them. Elements
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
values at the same points in both

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

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 Source

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($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

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download