【问题标题】:Extract shapefile value to point with R提取 shapefile 值以指向 R
【发布时间】:2023-03-26 12:50:02
【问题描述】:

我想在特定位置提取 shapefile 值。我使用的 shapefile 可以在here 找到,点击Download IHO Sea Areas 下载。形状文件包含所有可能的海洋。

我可以阅读并使用:

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

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

但是,我想提取特定位置的海值,比如说

p <- c(-20, 40)

【问题讨论】:

  • 那时您期望什么样的“价值”?海拔?但根据定义,这应该为零。还有其他价值吗?

标签: r gis rgdal


【解决方案1】:

可能有一种更简单的方法,但这是我的看法

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

输出点 x=-20, y=40

   xx                 NAME ID Gazetteer_ id
35  1 North Atlantic Ocean 23       1912 35

【讨论】:

    【解决方案2】:

    您可以使用 sp 包中的 over 函数:

    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)
    

    甚至更短:

    pts %over% shape
    

    【讨论】:

    • 这是最好的答案,100%。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2012-08-31
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多