rosscova - 14 days ago 6x
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.

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)
``````

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

Source (Stackoverflow)