从经度/纬度网格数据中提取 shapefile答案

【问题标题】：Extract shapefiles from longitude/latitude gridded data从经度/纬度网格数据中提取 shapefile
【发布时间】：2017-02-03 09:52:50
【问题描述】：

我有一些地中海海面温度值的网格数据，我已对其应用了聚类。我有 420 个具有三列结构（long、lat、value）的文件。特定文件的数据如下图所示

现在我想将聚类区域提取为 shapefile 以进行后处理。我找到了这篇文章 (https://gis.stackexchange.com/a/187800/9227) 并尝试像这样使用它的代码

# Packages
library(sp)
library(rgdal)
library(raster)

# Paths
ruta_datos<-"/home/meteo/PROJECTES/VERSUS/OUTPUT/DATA/CLUSTER_MED/"

setwd("~/PROJECTES/VERSUS/temp")

# File list
files <- list.files(path = ruta_datos, pattern = "SST-cluster-mitja-mensual")

for (i in 1:length(files)){

  datos<-read.csv(paste0(ruta_datos,files[i],sep=""),header=TRUE)
  nclusters<-max(datos$cluster)

  for (j in 1:nclusters){
    clust.dat<-subset(datos, cluster == j)

    coordinates(clust.dat)=~longitud+latitud

    proj4string(clust.dat)=CRS("+init=epsg:4326")
    pts = spTransform(clust.dat,CRS("+init=epsg:4326"))

    gridded(pts) = TRUE

    r = raster(pts)
    projection(r) = CRS("+init=epsg:4326")

    # make all values the same. Either do
    s <- r > -Inf

    # convert to polygons
    pp <- rasterToPolygons(s, dissolve=TRUE)

    # save shapefile
    shname<-paste("SST-shape-",substr(files[i],27,32),"-",j,sep="")
    writeOGR(pp, dsn = '.', layer = shname, driver = "ESRI Shapefile")    

  }

}

但代码停止并显示此错误消息

网格化(pts) = TRUE 建议的最小公差：1
points2grid 中的错误（点，公差，圆形）：尺寸 2 : 坐标间隔不是恒定的警告消息：在 points2grid(points, tolerance, round) 中：网格为空维度 1 中的列/行

我不明白在某个文件中它说坐标间隔不是恒定的，而它们确实是恒定的，从中派生聚类的原始 SST 数据位于全球的规则网格上。所有集群数据文件的大小相同，均为 4248 点。样本数据文件可用here

容差建议是什么意思？我一直在寻找解决方案，发现了一些使用SpatialPixelsDataFrame 的建议，但找不到如何申请。

任何帮助将不胜感激。谢谢。

【问题讨论】：

我尝试调试您的代码，并注意到j = 4 与您的示例数据发生错误。看看是不是在gridded(pts) = TRUE前加if (j == 4) debugonce(points2grid)，然后单步执行代码，能不能帮你找到问题？在 RStudio IDE 中效果最佳

标签： r shapefile r-raster sp

【解决方案1】：

我不是地理空间数据专家，但对我而言，如果您按集群过滤，数据确实不在网格上。据我了解，你从一个网格开始（凸集定期遥远的点）。

我尝试对您的代码进行以下修改并生成了一些文件，但我无法测试它们是否正确。原则是在所有数据上建立网格，然后在调用raster之前只对集群进行过滤。

这给出了：

files <- list.files(path = ruta_datos, pattern = "SST-cluster-mitja-mensual")

for (i in 1:length(files)){

  datos<-read.csv(paste0(ruta_datos,files[i],sep=""),header=TRUE)
  nclusters<-max(datos$cluster)

  for (j in 1:nclusters){
## clust.dat<-subset(datos, cluster == j)
clust.dat <- datos

coordinates(clust.dat)=~longitud+latitud

proj4string(clust.dat)=CRS("+init=epsg:4326")
pts = spTransform(clust.dat,CRS("+init=epsg:4326"))

gridded(pts) = TRUE

## r = raster(pts)
r= raster(pts[pts$cluster==j,])

projection(r) = CRS("+init=epsg:4326")

# make all values the same. Either do
s <- r > -Inf

# convert to polygons
pp <- rasterToPolygons(s, dissolve=TRUE)

# save shapefile
shname<-paste("SST-shape-",substr(files[i],27,32),"-",j,sep="")
writeOGR(pp, dsn = '.', layer = shname, driver = "ESRI Shapefile")    
  }
}

所以，注释中的两行并仅更新下面的行。

【讨论】：