rosscova - 8 months ago 39

R Question

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`

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

`[,,1]`

`[,,2]`

`[,,3]`

`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,]`

`> map[1,1,]`

[1] 1 50 100

I assume the information I need is contained in the

`mapObject`

I've been trying to figure out what the values in

`mapObject`

`[[1]]`

`[[2]]`

`x`

`y`

`NA`

`x`

`y`

Can anyone help me decipher the meaning of the contents of

`mapObject`

The

`map`

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

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

Source (Stackoverflow)