R - пакет «растровый» - RasterLayer («растр») не извлекается в SpatialPolygonsDataFrame («sp»); Просто возвращает объект NULL

Question

Используя пакет raster, у меня есть два набора данных, один растровый ASCII и один шейп-файл ESRI. Я хочу извлечь данные растра (температуры воды) в полную длину шейп-файла, который является береговой линией озера.

Форм-файл ESRI обрабатывается как SpatialPolygonsDataFrame при чтении с использованием функции shapefile().

shapefile <- shapefile("shore.shp",verbose=TRUE)

Я использовал функцию raster() для чтения в ASCII растра.

raster <- raster("1995_001.asc",native=TRUE,crs="+proj=merc +datum=WGS84 +ellps=WGS84 +units=m")

Координата справочной информации для шейпа является:

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

То растра (т.е. принуждаемой к следующему с помощью crs аргумента в функции raster()):

+proj=merc +datum=WGS84 +ellps=WGS84 +units=m +towgs84=0,0,0

Затем я использовал функцию spTransform() в пакете rgdal, чтобы обеспечить пространственную привязку шейп-файла к пространству растра.

spTransform(shapefile, CRS(projection(raster)))

Наконец, я представил следующее:

extract(raster,shapefile,method="simple",fun=mean,small=TRUE,na.rm=TRUE,df=FALSE)

Однако extract() просто возвращает NULL объект list типа. Я предполагаю, что этот вопрос подтверждается явным принуждением ссылок на координаты.

Кроме того, здесь приведены результаты использования функции show() на каждом наборе:

> show(raster) class : RasterLayer dimensions : 1024, 1024, 1048576 (nrow, ncol, ncell) resolution : 1800, 1800 (x, y) extent : -10288022, -8444822, 4675974, 6519174 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : layer values : -9999, 8.97 (min, max)

> show(shapefile) class : SpatialPolygonsDataFrame features : 1 extent : 597568.5, 998261.6, 278635.3, 668182.2 (xmin, xmax, ymin, ymax) coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 variables : 3 names : AREA, PERIMETER, HECTARES min values : 59682523455.695, 5543510.075, 5968252.346 max values : 59682523455.695, 5543510.075, 5968252.346

Я искал множество подобных вопросов на этих форумах, не имеющих разрешения. Может ли кто-нибудь дать мне (виртуальную) руку?

спасибо, что заранее.

Answer 1

Это (данные, которые не отображаются в пространственном перекрытии) является общей проблемой, вызванной путаницей относительно систем координат. Это может произойти по нескольким причинам.

1) регионы фактически не перекрывают друг друга

2) Возможно, ваше назначение АСБ к RasterLayer является неправильным (или что многоугольников). Пожалуйста, покажите объект (show(raster)); это почти всегда полезно при задании вопроса о пакете raster.

3) Вы не назначаете результат spTransform.Помните, что R обычно не выполняет модификации «на месте». Обязательно сделайте shapefile <- spTransform(shapefile, crs(raster))

Всегда выполняйте проверку на работоспособность и прокладывайте себе путь до тех пор, пока все не сработает. Место, чтобы начать здесь будет делать

plot(raster) 
plot(shapefile, add=TRUE) 

Чтобы увидеть, если какой-либо из полигонов участок поверх растровых данных.

Если это не сработает, посмотрите на экстенты. Например, как это (кстати, это также показывает, как можно создать автономный воспроизводимый пример/вопрос, что другие могут на самом деле ответ):

library(raster) 
# the spatial parameters of your raster 
r <- raster(ncol=100, nrow=100, xmn=-10288022, xmx=-8444822, ymn=4675974, ymx=6519174, crs="+proj=merc +datum=WGS84") 
values(r) <- 1:ncell(r) 

# the extent of your SpatialPolygons 
ep <- extent(597568.5, 998261.6, 278635.3, 668182.2) 
p <- as(ep, 'SpatialPolygons') 
crs(p) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83" 

# tranform both data set to longlat 
rgeo <- projectRaster(r, crs='+proj=longlat +datum=WGS84') 
pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84')) 

# plot on top of Google map 
library(dismo) 
e <- union(extent(rgeo), extent(pgeo)) 
g <- gmap(e,lonlat=TRUE) 
plot(g, inter=TRUE) 
plot(extent(rgeo), add=TRUE, col='red', lwd=2) 
plot(pgeo, add=TRUE, col='blue', lwd=2) 

Очевидно, что два источника данных делают не перекрываются. Только вы знаете, какая из двух может быть правильной, и теперь вы можете начать изучать, почему другой находится в неправильном месте.