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