Как рассчитать полигональные средства и отобразить их?

Question

У меня есть набор данных (pts), подобный этому:

x <- seq(-124.25,length=115,by=0.5)    
y <- seq(26.25,length=46,by=0.5)
z = 1:5290

longlat <- expand.grid(x = x, y = y)  # Create an X,Y grid
pts=data.frame(longlat,z) 
names(pts) <- c( "x","y","data")

Я знал, что могу отобразить фрейм данных (pts) на карту, выполнив:

library(sp)
library(rgdal)
library(raster)
library(maps)
coordinates(pts)=~x+y

proj4string(pts)=CRS("+init=epsg:4326") # set it to long, lat

pts = spTransform(pts,CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
pts <- as(pts, "SpatialPixelsDataFrame")
r = raster(pts)
projection(r) = CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

plot(r)
map("usa",add=T)

мой вопросэто то, как я могу рассчитать средние значения для 10 регионов EPA и отобразить средние значения?

регионы EPA можно найти в нижней части веб-страницы по адресу http://www.epa.gov/wed/pages/ecoregions/level_iii_iv.htm

Спасибо

Answer 1

Сначала прочитайте файл формы

er <- readOGR("Eco_Level_III_US.shp", "Eco_Level_III_US")

Затем убедитесь, что все растры r и экорегионы er имеют одинаковую проекцию

er.4326 <- spTransform(er, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

Извлеките растровые данные изшейп-файл (это может занять несколько минут), вычислите среднее значение и добавьте его к своим многоугольникам.

er.v <- extract(r, er.4326)
means <- sapply(er.v, mean)
er.4326$means <- means

И, наконец, построите его

 spplot(er.4326, "means")