Преобразование шейп-файла в растр

У меня возникла проблема с растрированием шейп-файла для создания точек на сетке 0,5 * 0,5. Шейп-файл представляет классификацию уровней риска (Низкий-0, Средний-100, Высокий-1000, Очень высокий-1500) глобальных коралловых рифов по отношению к комплексным угрозам.

Я взял код из другого примера, который отлично работает, но когда я пробую его для своих данных, я ничего не получаю от функции построения графика. См. ниже ссылку на шейп-файл и мой код:

Рифы под угрозой: глобальные комплексные угрозы

# Read shapefile into R
library(rgdal)
library(raster)    

int.threat.2030 <- readOGR(dsn = "Global_Threats/Integrated_Future", 
                           layer = "rf_int_2030_poly")

## Set up a raster "template" for a 0.5 degree grid
ext <- extent(-110, -50, 0, 35)
gridsize <- 0.5
r <- raster(ext, res=gridsize)

## Rasterize the shapefile
rr <- rasterize(int.threat.2030, r)

## Plot raster
plot(rr)

Есть идеи, где я могу ошибаться? Это проблема с самим шейп-файлом?

Пожалуйста и спасибо!


r shapefile rasterize
person RMFish    schedule 29.01.2016    source источник

Ответы (1)


arrow_upward
7
arrow_downward

Вы предположили, что многоугольники были в координатах долгота/широта, но это не так:

library(raster)
library(rgdal)
p <- shapefile('Global_Threats/Integrated_Future/rf_int_2030_poly.shp')
p

#class       : SpatialPolygonsDataFrame 
#features    : 63628 
#extent      : -18663508, 14601492, -3365385, 3410115  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=cea +lon_0=-160 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 3
#names       :    ID, THREAT, THREAT_TXT 
#min values  :     1,      0,   Critical 
#max values  : 63628,   2000,  Very High

Вы можете изменить проекцию

pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

а затем сделайте что-то вроде:

ext <- floor(extent(pgeo))
rr <- raster(ext, res=0.5)
rr <- rasterize(pgeo, rr, field=1)

Или сохраните исходную CRS и сделайте что-то вроде:

ext <- extent(p)
r <- raster(ext, res=50000)  
r <- rasterize(p, r, field=1)
plot(r)

Обратите внимание, что вы растрируете очень маленькие полигоны в большие растровые ячейки. Многоугольник считается «внутренним», если он покрывает центр ячейки (т. е. предполагается случай, когда многоугольники покрывают несколько ячеек). Таким образом, для этих данных вам нужно будет использовать гораздо более высокое разрешение (а затем, возможно, агрегировать результаты). В качестве альтернативы вы можете растеризовать центроиды полигонов.

Но ничего из вышеперечисленного на самом деле не имеет значения, поскольку вы делаете все это задом наперед. Полигоны явно получены из растра (посмотрите, какие они блочные), и растр доступен в наборе данных, на который вы указываете!

Итак, вместо растеризации выполните:

x <- raster('Global_Threats/Integrated_Future/rf_int_2030')
x
#class       : RasterLayer 
#dimensions  : 25456, 80150, 2040298400  (nrow, ncol, ncell)
#resolution  : 500, 500  (x, y)
#extent      : -20037508, 20037492, -6363885, 6364115  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 
#data source : C:\temp\Global_Threats\Integrated_Future\rf_int_2030 
#names       : rf_int_2030 
#values      : 0, 2000  (min, max)
#attributes  :
#   ID  COUNT THREAT_TXT
#    0  80971        Low
#  100 343535     Medium
# 1000 322231       High
# 1500 168518  Very High
# 2000  83598   Critical

Здесь изображена часть Палавана:

e <- extent(c(-8990636, -8929268, 1182946, 1256938))
plot(x, ext=e)
plot(p, add=TRUE)

Если вам нужно более низкое разрешение, см. raster::aggregate. Для другой системы отсчета координат см. raster::projectRaster.

person Robert Hijmans    schedule 30.01.2016
comment
RobertH Я не успел поблагодарить вас за помощь. Мой исходный код был основан на первом примере, который я встретил здесь в stackoverflow: stackoverflow.com/questions/14992197/ Разрешение данных в моем шейп-файле сделало невозможным растрирование до указанного мной разрешения. Мне очень понравилось ваше упрощение в конце поста. Ваше здоровье. - person RMFish; 17.03.2016
comment
Чтобы немного автоматизировать процесс, я создал github.com/brry/misc/blob/ master/shp2raster.R - person Berry Boessenkool; 25.11.2016