2017-02-03 20 views
1

У меня есть некоторые данные о значениях температуры поверхности моря в Средиземном море, к которым я применил кластеризацию. У меня 420 файлов с тремя структурами столбцов (long, lat, value). Данные для конкретного файла выглядит эта картаИзвлечь шейп-файлы из данных с координатами долготы/широты

enter image description here

Теперь я хочу, чтобы извлечь кластера области, как шейп для постобработки. Я нашел этот пост (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) = ИСТИНА предложил минимум допуска: 1
об ошибках в points2grid (точки, толерантность, круглые): измерение 2 : интервалы координаты не являются постоянными сообщениями Внимания: в points2grid (точки, толерантность, круглый): сетка имеет пустые столбцов/строки в размерности 1

Я не понимаю, что в определенном файле говорится, что интервалы координат не являются постоянными, пока они действительно есть, исходные данные SST, из которых была произведена кластеризация, находятся на регулярной сетке по всему миру. Все файлы данных кластера имеют одинаковый размер, 4248 точек. Файл пробных данных доступен here

Что означает предложение толерантности? Я искал решение и нашел предложение использовать SpatialPixelsDataFrame, но не смог узнать, как его применять.

Любая помощь будет оценена по достоинству. Благодарю.

+1

Я попытался отладить ваш код и заметил, что ошибка произошла для 'j' =' 4' с вашими данными примера. Посмотрите, добавляет ли 'if (j == 4) debugonce (points2grid)' before 'gridded (pts) = TRUE', а затем, перейдя через код, поможет вам найти проблему? Лучше всего работает в RStudio IDE –

ответ

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

Таким образом, две строки комментария и обновление только ниже линии.