1

我已经在 R-GIS 中工作了一段时间，现在我开始需要一个函数来计算像素面积（以公顷为单位），但栅格采用纬度/经度坐标和度数分辨率。

它有效，但如果有人可以改进某些部分以更普遍地使用，我将不胜感激。

# function to calculate the area in hectares
measureAreahec <- function(imageRaster) {
  R <- 6378.137                                # radius of earth in Km
  areaImagetemp <- area(imageRaster)
  xresVal <- xres(areaImagetemp)
  yresVal <- yres(areaImagetemp)
  coords <- data.frame(coordinates(areaImagetemp))
  colnames(coords) <- c("lon","lat")
  lon1 <- coords$lon - (xresVal/2)
  lon2 <- coords$lon + (xresVal/2)
  lat1 <- coords$lat - (yresVal/2)
  lat2 <- coords$lat + (yresVal/2)

  # width
  dLat <- (lat1-lat1)*pi/180
  dLon <- (lon2-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlongitude <- R * c * 1000

  # heigth
  dLat <- (lat1-lat2)*pi/180
  dLon <- (lon1-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlatitude <- R * c * 1000
  return (dlatitude * dlongitude / 10000)     #distance in meters}

# Applying the function
# raster 
r <- raster(nrow=180, ncol=360) 
projection(r)<- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

r <- setValues(r, rnorm(ncell(r), 1, 10))

vectorArea <- measureAreahec(r)

rasterArea <- setValues(r, vectorArea)
### end
4

0 回答 0