在 R 中计算高纬度栅格面积的最准确方法？

编程问答 2022-05-26

问题描述

我希望计算跨纬度的多个栅格（~110 公里乘 20 公里）的面积（以 m^2 为单位）。我无法就最佳方式达成一致，以处理经度与纬度之间关系的变化。我发现的三种不同方法包括将 raster::area 函数应用于未投影的栅格以计算平均像元面积，另一种将 raster::area 函数应用于从具有 +proj=eqearth 投影的栅格导出的多边形，最后将 rgeos:::gArea 函数应用于从具有 +proj=eqearth 投影的栅格导出的多边形。我没有得到类似的结果。我在下面提供了一些代码。最好的技术是什么？提前致谢！

library(raster)
library(rgeos)
library(sp)
library(rgdal)

#generate a polygon to start off

#coordinates for outside of donut polygon
outer <- polygon(rbind(c(-15.5376528679349,79.4881322107174),c(-7.48352810724904,79.908837268938),c(-6.66906605279766,78.6475974600555),c(-15.4471570841069,78.7892258712421),c(-15.5376528679349,79.4881322107174)))

#coordinates for inside of donut polygon
inner <- polygon(rbind(c(-13.6190362010361,79.3320228496385),c(-13.7164252183978,79.039680481742),c(-9.04175238503622,c(-9.23653041975962,79.3859727584367),c(-13.6190362010361,79.3320228496385)),hole = T)

#create polygon from coordinates
h1 <- polygons(list(outer,inner),"donut1")

# create spatial polygons object
onedonut <- Spatialpolygons(list(h1))

#equal area projection
projection(onedonut) <- CRS("+proj=eqearth")



Convert polygon to raster

r <- raster(ncol=100,nrow=100)
extent(r) <- extent(onedonut)
onedonut.raster <- rasterize(onedonut,r)


#Calculate area using raster::area on non-projected raster

   #raster area calculation
      #get sizes of all cells in raster
    cell_size_raster<-area(onedonut.raster,na.rm=TRUE,weights=FALSE)

    #delete NAs from vector of all raster cells
    cell_size_raster<-cell_size_raster[!is.na(onedonut.raster)]

    #compute area of all cells
    segment_area_raster <-length(cell_size_raster)*median(cell_size_raster)

#Calculating area using raster::area on projected shapefile

    cell_size_sp<-area(onedonut,weights=FALSE)

    #compute area [km2] of all polygons in shape file
    segment_area_equalarea <-length(cell_size_sp)*median(cell_size_sp)


#Calculate area using rgeos::gArea on projected shapefile

rgeos_shapefile_area <- gArea(onedonut)
  
  


#Are these values close?

segment_area_raster
segment_area_equalarea
rgeos_shapefile_area

#The values are not the same

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

projectionrrrasterrgeo-shapefilespatial