使用西南角将 xy 栅格化为栅格

Question

  library(terra)
  
  dat <- structure(list(lon = c(73.0754166666667, 73.0754166666667, 73.0754166666667, 
                                73.0754166666667, 73.07625, 73.07625, 73.0770833333333, 
                                73.0770833333333, 73.0770833333333, 73.0770833333333), 
                  lat = c(-0.590416666666666, -0.587916666666666, -0.587083333333332, -0.58625, 
                          -0.592083333333333, -0.589583333333333, -0.593749999999999, 
                          -0.592083333333333, -0.589583333333333, -0.588749999999999), 
                  value = c(0.03, 0.02, 0.06, 0.032, 0.043, 0.023, 0.0123, 0.051, 0.0441, 0.015)), 
                  row.names = c(NA, -10L), class = "data.frame")

我有一个名为 dat 的数据框，其中包含 lon 和 lat 以及一些与之相关的值。然而，经度和纬度 代表西南角而不是每行的质心。

我想要做的是创建一个栅格并用

dat

中的值填充它。为了做到这一点，我需要首先创建一个空栅格，其参数定义为。

  pixelsize <- 0.0008333333
  nr <- 6000
  nc <- 6000
  west_bound <- 70 # west_bound  and southbound represents southwest corner of the empty raster
  south_bound <- -5 

  # Calculate the extent of the raster
  east_bound <- west_bound + ncols * resolution
  north_bound <- south_bound + nrows * resolution
  
  blank_canvas <-
    rast(nrows = nr, #create a blank raster
         ncols = nc, 
         xmin = west_bound, 
         xmax = east_bound, 
         ymin = south_bound, 
         ymax = north_bound, 
         crs = 'epsg:4326',
         resolution = pixelsize)
  
  values(blank_canvas) <- -999
  
  dat_pts <- terra::vect(dat, geom = c("lon", "lat"), crs = 'epsg:4326') # convert to points
  r <- terra::rasterize(x = dat_pts, y = blank_canvas, field = "value", touches = T) 

但是，生成的光栅将

dat

中的每个点视为质心，而不是单个像元的西南角

    r_poly <- as.polygons(r)
    plot(r_poly)
    plot(dat_pts, add = T)

Answer 1

如果栅格和点匹配，则显然坐标指的是像元的中心。也许您没有正确定义光栅？

当然，如果你认为合理的话，你可以移动这些点。

library(terra)
d <- data.frame(lon = c(73.0754166666667, 73.0754166666667, 73.0754166666667, 
                     73.0754166666667, 73.07625, 73.07625, 73.0770833333333, 
                     73.0770833333333, 73.0770833333333, 73.0770833333333), 
                lat = c(-0.590416666666666, -0.587916666666666, -0.587083333333332, -0.58625, 
                     -0.592083333333333, -0.589583333333333, -0.593749999999999, 
                     -0.592083333333333, -0.589583333333333, -0.588749999999999), 
                value = c(0.03, 0.02, 0.06, 0.032, 0.043, 0.023, 0.0123, 0.051, 0.0441, 0.015))
        
x <- rast(ncols=6000, nrows=6000, xmin=70, xmax=75, ymin=-5, ymax=0, crs='EPSG:4326')

r <- res(x)
d$lon <- d$lon + r[1]/2
d$lat <- d$lat + r[2]/2

xd <- rasterize(d[,1:2], x, d[,3])