CODEWITHSUNDEEP
rshapefiler-rasterr-sp
pacomet
pacomet

Reputation: 5141

Extract shapefiles from longitude/latitude gridded data

I have some gridded data of sea surface temperature values in the Mediterranean to which I've applied clustering. I have 420 files with three columns structure (long,lat,value). The data for a particular file looks like this map

enter image description here

Now I want to extract the cluster areas as shapefile for postprocessing. I have found this post (https://gis.stackexchange.com/a/187800/9227) and tried to use its code like this 

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

  }

}

But the code stops for with this error message

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

I don't understand that at a certain file it says that coordinate intervals are not constant while they indeed are, original SST data from which clustering was derived are on a regular grid over the whole globe. All cluster data files have the same size, 4248 points. A sample data file is available here

What does the tolerance suggestion means? I've been looking for a solution and found some suggestion to use SpatialPixelsDataFrame but couldn't find out how to apply.

Any help would be appreciated. Thanks.

Upvotes: 1

Views: 511

Answers (1)

Eric Lecoutre
Eric Lecoutre

Reputation: 1481

I am not an expert of geospatial data but for me, if you filter on cluster, data are indeed not on a grid. So far as I understand, you start from a grid (convex set of regularly distant points).

I tried following modifications to your code and some files are generated but I can't test whether they are correct or not. Principle is to build the grid on all data then only filter on cluster before calling raster.

This gives:

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

So, two lines in comment and update just the line below.

Upvotes: 1

Related Questions