преобразовать шейп-файл нескольких полигонов в растр в R - PullRequest
Новая
       83

преобразовать шейп-файл нескольких полигонов в растр в R

0 голосов
/ 02 июня 2018

У меня большие проблемы с преобразованием полигонов в растр в R. Я хотел сделать следующее: у меня есть шейп-файл (то есть полигоны) из 574 видов.То есть в таблице атрибутов она имеет 574 строки (т. Е. FID находится в диапазоне от 0 до 573).Подмножество данных можно найти здесь: https://drive.google.com/file/d/1AdTChjerCXopE1-PZIAPp5ZXISdq45i8/view?usp=sharing

Я хотел преобразовать его в растр.В выходном растре я вижу, что минимальное и максимальное значения равны 1 и 574. Я подозреваю, что это получение идентификатора поля в ячейке в виде значения пикселя, что не должно.Значения ячеек должны исходить от покрытия полигонов.Любая помощь будет высоко оценена.Ниже приведен пример кода: 

library(raster)
library(rgdal)
library(maptools)


# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

# load species polygons for all 574 species
# already in the home directory

sp <- readShapePoly("AllP574sp.shp", proj4string = projection1)

# load a raster file to use as a mirror for rasterize
raster1 <- raster("/data/projects/MeanTemp2050rcp45_BCC_CSM1_1.tif")


r.sp <- rasterize(sp, raster1) # rasterize our species polygon to the same resoluton of loaded raster

writeRaster(r.sp, "/data/projects/all574spRaster", format = "GTiff", overwrite = TRUE)

Свойства выходного растра приведены ниже: 

class       : RasterLayer 
dimensions  : 18000, 43200, 777600000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
data source : E:\all574spRaster.tif 
names       : all574spRaster 
values      : 1, 574  (min, max)

1 Ответ

0 голосов
/ 02 июня 2018

Когда вы использовали функцию rasterize, важно указать аргумент field, иначе по умолчанию он попытается создать его для вас;похоже, что в вашем случае он создал один из столбца FID.

Я сделал несколько предположений, чтобы восстановить рабочий набор полигонов, который может быть похож на ваш. 

library(maptools)
library(rgdal)
library(sp)
library(geosphere)

# set seed for duplicatable results
set.seed(1)

# some data that looks a little like yours
BINOMIAL <- c("Controversial chimneyswift", "Dull dungbeetle", 
"Easternmost eel", "Jumping jaeger", "Qualified queenconch")
FID <- 0:(length(BINOMIAL) - 1)
RANGE <-  runif(length(BINOMIAL), min = 118, max = 3875370)
MyData <- cbind.data.frame(FID, BINOMIAL, RANGE)
row.names(MyData) <- FID

# some semi-random polygons in your extent box
ext <- extent(c(-180, 180, -60, 90))
create_polygon <- function(n = 4, lat, lon, r) {
  lengths <- rnorm(n, r, r/3)
  smoother_lengths <- c(sort(lengths), rev(sort(lengths)))
  lengths <- smoother_lengths[sort(sample(n * 2, n))]
  lengths <- rep(lengths[1], length(lengths))
  directions <- sort(runif(n, 0, 360))
  p <- cbind(lon, lat) 
  vertices <- t(mapply(destPoint, b = directions, 
                       d = lengths, MoreArgs = list(p = p)))
  vertices <- rbind(vertices, vertices[1, ])
  sapply(vertices[,1], min, ext@xmax)
  sapply(vertices[,1], max, ext@xmin)
  sapply(vertices[,2], min, ext@ymax)
  sapply(vertices[,2], max, ext@ymin)
  Polygon(vertices) 
}
rand_lats <- runif(nrow(MyData), min = -50, max = 60)
rand_lons <- runif(nrow(MyData), min = -100, max = 100)
rand_sides <- sample(4:20, nrow(MyData), replace = TRUE)
rand_sizes <- rnorm(nrow(MyData), mean = 5e+06, sd = 1e+06)
make_species_polygon <- function(i) {
  p.i <- list(create_polygon(rand_sides[i], rand_lats[i], 
                             rand_lons[i], rand_sizes[i]))
  P.i <- Polygons(p.i, FID[i])
}

polys <- SpatialPolygons(lapply(1:nrow(MyData), make_species_polygon))
spdf <- SpatialPolygonsDataFrame(Sr = polys, data = MyData)

t.shp <- tempfile(pattern = "MyShapefile", fileext = ".shp")
raster::shapefile(spdf, t.shp)

На данный моментв вашем временном каталоге есть шейп-файл, имя которого хранится в переменной t.shp.Я намерен, чтобы этот шейп-файл был работоспособной копией любого большого шейп-файла, который является вашим истинным.Итак, теперь мы можем посмотреть на ваш код, что он делал и что бы вы хотели, чтобы он делал: 

## now we get into your code
library(raster)
library(rgdal)
library(maptools)

# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

## 
## I don't know what your shapefile looks like exactly, 
## but substituting `t.shp` the tempfile that I created above
## also since the function readShapePoly is deprecated
## instead I use the recommended new function rgdal::readOGR() 
## 
sp <- rgdal::readOGR(t.shp)

##
## I don't know what your tiff file looks like exactly,
## but I can duplicate its characteristics
## for speed I have decreased resolution by a factor of 10
##
raster1 <- raster(nrow = 1800, ncol = 4320, ext)

# rasterize our species polygon to the same resoluton of loaded raster
r.sp <- rasterize(x = sp, y = raster1, field = MyData$RANGE) 

t.tif <- tempfile(pattern = "MyRastfile", fileext = ".tif")
writeRaster(r.sp, t.tif, format = "GTiff", overwrite = TRUE)

со следующим результатом: 

raster(t.tif)
class       : RasterLayer 
dimensions  : 1800, 4320, 7776000  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : [["a file name in your temp directory"]] 
names       : MyRastfile1034368f3cec 
values      : 781686.3, 3519652  (min, max)

Результаттеперь отображаются значения, которые взяты из столбца RANGE вместо столбца FID.

...