6

I want to extract a shapefile value at a particular location. The shapefile I am using is can be found here, downloaded by clicking Download IHO Sea Areas. The shape file contains all the possible seas.

I can read it and plot it using:

require("maptools")
require(rgdal)
require(sp)

ogrListLayers("World_Seas.shp")
shape <- readOGR("World_Seas.shp", layer="World_Seas")

However, I would like to extract the sea value for a particular location, say

p <- c(-20, 40)
metasequoia
  • 7,014
  • 5
  • 41
  • 54
user3910073
  • 511
  • 1
  • 6
  • 23
  • What kind of 'value' do you expect at that point? Elevation? But that should be zero by definition. Any other value? – Felix Dec 14 '15 at 20:32

2 Answers2

4

there is a probably an easier way but here is my take

require("maptools")
require(rgdal)
require(sp)
library(plyr)
library(dplyr)

setwd("/Users/drisk/Downloads/seas")
ogrListLayers("World_Seas.shp")
shape=readOGR("World_Seas.shp", layer="World_Seas") 

datapol <- data.frame(shape)
pointtoplot <- data.frame(x=-20, y=40)
coordinates(pointtoplot) <- ~ x + y 
proj4string(pointtoplot) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#
#function over from package sp
test <- data.frame(xx=over(shape, pointtoplot))
combine <- cbind(test, datapol)
combine <- na.omit(combine) #only one point left

output for your point x=-20, y=40

   xx                 NAME ID Gazetteer_ id
35  1 North Atlantic Ocean 23       1912 35
MLavoie
  • 9,671
  • 41
  • 36
  • 56
2

You could use the over function from the sp package:

library(rgdal)
library(sp)
library(raster)

shape <- shapefile("~/tmp/World_Seas.shp")
head(shape)

plot(shape[shape$ID == 35, ], axes = TRUE)
points(pts)

pts <- SpatialPoints(cbind(-20, 40), 
                     proj4string = CRS(proj4string(shape)))
over(pts, shape)

or even shorter:

pts %over% shape
johannes
  • 14,043
  • 5
  • 40
  • 51