Извлеките Растр из Растра, используя шейп-файл Polygon в R

13

Я новичок в R и использую растровый пакет. У меня проблема с извлечением полигонов из существующего растрового файла. Если я использую

extract(raster, poly_shape)

Функция на растре всегда создает список с данными. Что я действительно хочу, так это извлечь еще один растровый файл, который я могу снова загрузить с ArcGIS. Прочитав немного больше, я думаю, что функция обрезки - это то, что мне действительно нужно. Но когда я пытаюсь использовать эту функцию

crop(raster, poly_shape)

Я получаю эту ошибку:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

Файлы raster и poly_shape одинаковы для обеих функций. Можете ли вы сказать мне, что здесь может быть не так? Верно ли, что функция кадрирования создает другой растр, а не список?

РЕДАКТИРОВАТЬ : функция экстента () не работает для меня. Я все еще получаю ту же ошибку. Но я уверен, что 2 набора данных перекрываются! С

extract(raster, poly_shape)

Я получаю правильные данные из него. Просто как список, а не как растр, как я хочу его иметь. Я только что загрузил наборы данных в ArcGIS раньше, и они очень хорошо подходят, поэтому я не проверял проекцию. Сейчас попробовал

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

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

  • Я вырезал точную часть многоугольника, которую я извлек в R также в ArcGIS
  • Я рассчитал сумму всех значений извлеченного R полигона (список)
  • Я рассчитал сумму всех растровых ячеек, которые я вырезал в ArcGIS

2 имеют точно такой же результат, поэтому я думаю, что следует сделать вывод, что функция извлечения работала правильно. Теперь у меня есть 2 варианта, я думаю:

  1. Мне нужен способ вытащить растр из извлеченного списка снова или
  2. 2 набора данных (растр + poly_shape) должны использовать одну и ту же проекцию, и функция обрезки должна работать

Что бы вы предложили сделать здесь?

Lars
источник
Что если это RGBI растр с 4 полосами? Полосы потеряны до сих пор ...
Дорис
Добро пожаловать в ГИС ЮВ! Как новый пользователь, пожалуйста, ознакомьтесь с коротким туром . Затем рассмотрите возможность редактирования своего ответа, чтобы предоставить дополнительную информацию и ссылки. Смотрите, как мне написать хороший ответ? для получения дополнительной информации.
Энди
Если у вас есть новый вопрос, задайте его, нажав кнопку « Задать вопрос» . Включите ссылку на этот вопрос, если это помогает обеспечить контекст. - Из обзора
Эрик

Ответы:

19

Функция извлечения ведет себя точно так, как и должна. Вы можете заставить функцию обрезки использовать экстент многоугольника, а затем замаскировать объект, чтобы получить точный растр, представляющий область многоугольника. Если вы продолжаете получать сообщение об ошибке, это означает, что ваши данные фактически не перекрываются.

Пожалуйста, имейте в виду, что R не выполняет проекцию "на лету", поэтому проверьте свои прогнозы. Вы можете проверить, перекрываются ли ваши экстенты, используя функцию "stretch ()".

Вот пример обрезки с использованием экстента многоугольника с последующим маскированием полученного растра с использованием «растеризованного» многоугольника.

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)
Джеффри Эванс
источник
4
extract () выполняет перепроецирование на лету, но crop () - нет. Это может объяснить некоторую путаницу
mdsumner
@Jefferey crop () и mask () обрезают растр только в соответствии с прямоугольными экстентами многоугольника, но не обрезают его из-за границы многоугольника. Любая идея, какие команды могут обрезать растр в границах данного полигона?
Чеф
1
@Chintan Sheth, для подмножества маски в пределах многоугольника необходимо иметь растр, представляющий значения внутри многоугольника. Вот почему вы растеризуете полигон подмножества, а затем маскируете его. Шаг кадрирования состоит в том, чтобы уменьшить экстент растра до того же размера, что и в подмножестве полигонов, что в прошлом, если пропустить его, приводило к ошибке несоответствия экстентов.
Джеффри Эванс
spTransformиз spпакета (который иногда автоматически загружается с другими пакетами пространственного R) позволяет перепроецировать так, чтобы оба файла были в одной проекции, например. good_poly=spTransform(spolygon, CRSobj=crs(raster_file))
user3386170
@ user3386170, а? Не уверен, к чему ты клонишь. Этот вопрос возник в то время, когда растровый пакет просто добавил «проекцию на лету» в некоторые из его функций. Эти функции ранее выдавали ошибку, когда проекции не совпадали, но это сообщение было от 2014 года. Вы также должны помнить, что всегда следует загружать rgdal при использовании sp :: spTransform (), так как он добавляет дополнительные важные функции за кулисы.
Джеффри Эванс
8

На самом деле я искал mask()функцию.

mask(raster, poly_shape)

работает без ошибок и возвращает то, что искал.

Lars
источник
2
Перепроектируйте ваши данные так, чтобы они находились в том же пространстве проекции. Даже в ArcGIS, где проекция на лету автоматическая, очень плохо проводить анализ в разных проекциях. С данными в разных проекциях ваши данные не будут иметь общих экстентов, что является ошибкой, которую вы получаете.
Джеффри Эванс
Используйте projectExtent (), чтобы получить только степень обрезки растра.
mdsumner
Чтобы соответствовать формату вопросов и ответов на сайте, его следует поместить в основную часть вопроса в виде правки / обновления (а затем добавить комментарий к ответу, на который есть «ответ», чтобы они знали, что есть еще на что посмотреть).
Мэтт Уилки
@mattwilkie Извините, что не соответствовал формату, но мой текст был слишком длинным, чтобы разместить его здесь в качестве комментария. @JeffreyEvans Я попробовал следующее: projection(raster) = projection(poly_shape)и наоборот , projection(poly_shape) = projection(raster)но оба способа производят ту же ошибку: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. Есть ли способ, как я могу увидеть, какая проекция используется на лету, используя функцию extract () (потому что эта, очевидно, работает)?
Lars
1
На самом деле я искал функцию mask (). mask(raster, poly_shape)работает без ошибок и возвращает то, что искал.
Ларс
-1

Экстент работает просто отлично ... Я думаю, что Xmin, Xmax, Ymin и Ymax вашего экстента отличаются от X и Y вашего растра - т.е. они установлены напротив

ToNoY
источник