Я пытаюсь извлечь пропорциональные значения земного покрова из нескольких шейп-файлов USFWS National Wetlands Inventory в буферных зонах вокруг фрейма данных пространственных точек.При построении в QGIS шейп-файлы содержат множество полигонов, представляющих водно-болотные угодья, каждый из которых имеет свою классификацию.Моя цель состоит в том, чтобы преобразовать этот шейп-файл в растровый слой формального класса и извлечь эти значения с помощью raster :: extract ().
Я рассчитал пропорциональное покрытие земли для буферизованных точек успешно в унисон с библиотекой (FedData) и get_nlcd ();Тем не менее, я был полностью неудачным при попытке работать с ESRI .shp файлами.Каждый раз, когда я пытался следовать чужим инструкциям по растеризации шейп-файлов, извлеченный фрейм данных полон нулевых значений.Я не уверен, пропустил ли я шаг, или это просто форма файла, с которым я работаю, но я надеялся, что кто-то может предложить решение.
library(sp)
library(raster)
library(rgdal)
#making a spatial data frame from coordinates
coords = cbind(SeniorThesisSamples$Longitude, SeniorThesisSamples$Latitude)
points_spdf = SpatialPoints(coords, proj4string = crs("+proj=aea +ellps=GRS80 +datum=NAD83"))
#Extending spatial extent to make sure all points can have an X km buffer -----------------
r = raster(ncol = 3, nrow = 3)
extent(r) = extent(points_spdf) * 1.5
crs(r) = crs(points_spdf)
#read in National Wetlands Inventory ESRI shapefile
RI.nwi <- raster::shapefile("~/Documents/Senior Thesis Stuff/NWI Shapefiles/RI_shapefile_wetlands/RI_Wetlands.shp")
# Transform .shp to match CRS of the spdf
RI.nwi <- spTransform(RI.nwi, crs(r))
#creating an empty raster with the extent of my spdf and setting resolution to 100m
empty.raster <- raster(extent(r), res = 100)
crs(empty.raster) = crs(r)
RI.nwi.raster <- rasterize(RI.nwi, empty.raster)
#extracting raster values at a set buffer distance and put into data frame
wetlands1km <- raster::extract(x = RI.nwi.raster, y = points_spdf, buffer = 1000, df = T)
#this part I know works because I have used it effectively to extract NLCD data
detach(package:plyr) # this is so the n() works
prop.wetlands1km <- as.data.frame(wetlands1km) %>%
setNames(c("ID", "class")) %>%
group_by(ID, class) %>%
summarise(n = n()) %>%
mutate(pland = n / sum(n)) %>%
ungroup() %>%
dplyr::select(ID, class, pland) %>% # keep only these vars
tidyr::complete(ID, nesting(class), fill = list(pland = 0)) %>% # Fill in implicit landcover 0s
tidyr::spread(class, pland)
Код выполняется:но у меня остался пустой фрейм данных, который должен (из предыдущего опыта работы с данными NLCD) иметь классификацию полей в качестве заголовков столбцов и процентных значений в каждой ячейке.
Я загрузил несколько шейп-файлов уровня состояния по этой ссылке (https://www.fws.gov/wetlands/Data/State-Downloads.html) для справки. Точки в моем spdf простираются за пределы только Род-Айленда, но это был простейший шейп-файл, потому что он очень маленький. Бонусные баллы могут помочь любому найти способ объединить несколько из этихсоставьте шейп-файлы вместе (ME, NH, MA, RI, CT, NY, NJ, DE, MD, VA), чтобы я мог запустить извлечение значений за один раз! Любая помощь очень ценится!
РЕДАКТИРОВАТЬ:Прикреплен пример запрашиваемого кода
> dput(head(SeniorThesisSamples))
structure(list(SiteName = c("Quonochontaug", "Quonochontaug",
"Quonochontaug", "Quonochontaug"), State = c("RI", "RI", "RI",
"RI"), SiteID = c("220450_p3", "220450_p3", "220450_p3", "220450_p3"
), Latitude = c(41.3335569999999, 41.3335569999999, 41.3335569999999,
41.3335569999999), Longitude = c(-71.7167649999999, -71.7167649999999,
-71.7167649999999, -71.7167649999999), Species = c("SALS", "SALS",
"SALS", "SALS"), BandNum = c(278184716, 278184718, 278184721,
278184715), HgConcentration = c(0.422, 0.32, 0.374, 0.309), Month = c(7,
7, 7, 7), Day = c(29, 29, 29, 29), Year = c(2018, 2018, 2018,
2018), Date = structure(c(17741, 17741, 17741, 17741), class = "Date")), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -4L))
> RI.nwi
class : SpatialPolygonsDataFrame
features : 55398
extent : 1960230, 2051242, 2271966, 2388062 (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 : 5
names : ATTRIBUTE, WETLAND_TY, ACRES, Shape_Leng, Shape_Area
min values : E1AB1L, Estuarine and Marine Deepwater, 1.000117e+00, 1.000005e+02, 1.000025e+03
max values : R5UBH, Riverine, 9.998890e-03, 9.999467e+01, 9.999969e+02
> head(RI.nwi)
ATTRIBUTE WETLAND_TY ACRES Shape_Leng Shape_Area
0 E2ABM Estuarine and Marine Wetland 0.5421139 192.1974 2193.857
1 E2ABM Estuarine and Marine Wetland 0.4598518 197.5223 1860.954
2 E2ABM Estuarine and Marine Wetland 2.9291088 453.6280 11853.683
3 E2AB3M Estuarine and Marine Wetland 2.0402979 514.5920 8256.793
4 E1AB1L Estuarine and Marine Deepwater 6.3219543 780.8088 25584.041
5 E1AB1L Estuarine and Marine Deepwater 2.0643299 361.9401 8354.047