从 longitude/latitude 网格化数据中提取 shapefile

Extract shapefiles from longitude/latitude gridded data

我有一些应用了聚类的地中海海面温度值的网格化数据。我有 420 个具有三列结构(长、纬度、值)的文件。特定文件的数据如下图所示

现在我想将聚类区域提取为 shapefile 以供 post 处理。我找到了这个 post (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")    

  }

}

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

gridded(pts) = TRUE suggested tolerance minimum: 1
Error in points2grid(points, tolerance, round) : dimension 2 : coordinate intervals are not constant Warning message: In points2grid(points, tolerance, round) : grid has empty column/rows in dimension 1

我不明白在某个文件中它说坐标间隔不是恒定的而它们确实是,从中导出聚类的原始SST数据在整个地球上处于规则网格上。所有簇数据文件的大小相同,均为 4248 点。样本数据文件可用 here

宽容建议是什么意思?我一直在寻找解决方案并找到了一些使用 SpatialPixelsDataFrame 的建议,但找不到如何应用。

如有任何帮助,我们将不胜感激。谢谢。

我不是地理空间数据专家,但对我来说,如果你在集群上过滤,数据确实不在网格上。据我所知,你从一个网格(一组规则远点的凸集)开始。

我尝试对您的代码进行以下修改并生成了一些文件,但我无法测试它们是否正确。 原则是在所有数据上构建网格,然后在调用 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")    
  }
}

因此,注释两行并仅更新下面的行。