我的 CRS 指定的地图单位是什么？
What is the map unit specified by my CRS?
我对 R 的空间工作还比较陌生，最近 运行 遇到了以下问题：
我有气候数据网格化（最高温度，每月），想提取某些位置的观测值，使用半径为 20KM 的圆形缓冲区。让我先设置并绘制它...
library(raster)
# load NetCDF file into SpatialRaster
ncin <- raster::brick(x='cru_ts4.05.1901.2020.tmx.dat.nc')
#> Warning in .varName(nc, varname, warn = warn): varname used is: tmx
#> If that is not correct, you can set it to one of: tmx, stn
# get shape file
poly <- readRDS("C:/Users/m/shapes/gadm36_NGA_1_sf.rds") %>%
st_transform(proj4string(ncin))
ncin <- raster::crop(ncin, poly, snap = "out")
#clip raster outside the poly border
# result is a raster brick with dimensions (lat=20, lon=25, t=1440)
load ("C:/Users/m/NGA_panel.Rda") # load household data, with coordinates
# plot households as dots on the grid; see whether all parts display correctly!
plot(ncin$X2020.12.16, main="HH-locations")
# plot(poly$geometry, add = TRUE)
points(NGA_panel$LON_DD_MOD, NGA_panel$LAT_DD_MOD, pch=19, cex=0.5, col =1)
Here 是结果图。黑点代表家庭。对于它们中的每一个，我想使用半径为 20 KM 的圆形缓冲区从网格中提取一个值。提取的值是缓冲区内栅格值的平均值。
我的 CRS 是 +proj=longlat +datum=WGS84 +no_defs，根据 documentation of raster::buffer()，它不一定将米作为其默认地图单位。 st_crs(ncin) 产生以下详细信息：
st_crs(ncin) # does the current crs take meters as length-unit?
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
以上 LENGTHUNIT["metre",1] 是否意味着缓冲区半径将以米为单位正确缩放？继续提取让我...
NGA_spatial <- SpatialPointsDataFrame(
NGA_panel[108:107], proj4string=ncin@crs, NGA_panel) # spatialize dataset (i.e. make SPDF)
NGA_spatial <- NGA_spatial[sample(nrow(NGA_spatial), 100), ] #### LIMIT to 100 random rows FOR TESTING!
# extract values for households, using a circular buffer with radius 20KM:
data <- raster::extract(ncin,
NGA_spatial, # SPDF with centroids for buffer
buffer = 20000, # buffer size, units depend on CRS (meters or degrees?)
fun=mean, # what value to extract
df=TRUE) # create a data frame?
summary(data[sample(ncol(data), 3)]) # summary statistics of 3 random columns
#> X2005.09.16 X2000.02.15 X2001.02.15
#> Min. :28.80 Min. :29.20 Min. :31.20
#> 1st Qu.:29.77 1st Qu.:30.30 1st Qu.:32.17
#> Median :30.70 Median :33.30 Median :34.10
#> Mean :31.21 Mean :32.46 Mean :33.63
#> 3rd Qu.:32.25 3rd Qu.:33.80 3rd Qu.:34.55
#> Max. :37.10 Max. :35.00 Max. :36.00
显然我提取的值存在一些变化，但我不确定这是通过 20000 米的缓冲区（如我所愿）还是 20000 度的缓冲区（这很荒谬）提取的。
有人能告诉我这是否正确地以米为单位提取，是否有办法控制它（例如通过某种方式绘制缓冲区）？非常感谢！
您需要查看WKT的分组：
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
里面有 3 个 *UNIT 块：
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]
和 2 个轴：
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
第一个适用于椭圆体（以米为单位），另外两个是您感兴趣的，因为它们定义了数据的存储方式（和呈现方式）以及以度为单位（与未投影的预期一样） WGS CRS）。因此，您的缓冲区将以度为单位，正如您所说的那样是荒谬的，只要没有任何中断，它将 return 您的所有数据。
要以米为单位工作，您需要将数据重新投影到平面本地 CRS（最好是对您的数据集等距离的 CRS）。
?raster::extract 参数“buffer”的文档指出
If the data are not projected (latitude/longitude), the unit should be
meters. Otherwise it should be in map-units (typically also meters).
+proj=longlat +datum=WGS84明明是(?)latitude/longitude，还是用米吧
如下所示，取自 ?extract 中的示例（略有改动）：
library(raster)
r <- raster(ncols=36, nrows=18, vals=1:(18*36))
xy <- cbind(0, seq(20, 80, by=20))
xy
# [,1] [,2]
#[1,] 0 20
#[2,] 0 40
#[3,] 0 60
#[4,] 0 80
extract(r, xy, buffer=1000000)
#[[1]]
#[1] 234 235 270 271
#[[2]]
#[1] 162 163 198 199
#[[3]]
#[1] 89 90 91 92 126 127
#[[4]]
# [1] 13 14 15 16 17 18 19 20 21 22 23 24 51 52 53 54 55 56 57 58
如您所见，当您向极点移动时，您会获得越来越多的值（即缓冲区内的单元格）。这是预期的，因为纵向距离减小，大致与 cos(latitude).
成正比
我对 R 的空间工作还比较陌生，最近 运行 遇到了以下问题：
我有气候数据网格化（最高温度，每月），想提取某些位置的观测值，使用半径为 20KM 的圆形缓冲区。让我先设置并绘制它...
library(raster)
# load NetCDF file into SpatialRaster
ncin <- raster::brick(x='cru_ts4.05.1901.2020.tmx.dat.nc')
#> Warning in .varName(nc, varname, warn = warn): varname used is: tmx
#> If that is not correct, you can set it to one of: tmx, stn
# get shape file
poly <- readRDS("C:/Users/m/shapes/gadm36_NGA_1_sf.rds") %>%
st_transform(proj4string(ncin))
ncin <- raster::crop(ncin, poly, snap = "out")
#clip raster outside the poly border
# result is a raster brick with dimensions (lat=20, lon=25, t=1440)
load ("C:/Users/m/NGA_panel.Rda") # load household data, with coordinates
# plot households as dots on the grid; see whether all parts display correctly!
plot(ncin$X2020.12.16, main="HH-locations")
# plot(poly$geometry, add = TRUE)
points(NGA_panel$LON_DD_MOD, NGA_panel$LAT_DD_MOD, pch=19, cex=0.5, col =1)
Here 是结果图。黑点代表家庭。对于它们中的每一个，我想使用半径为 20 KM 的圆形缓冲区从网格中提取一个值。提取的值是缓冲区内栅格值的平均值。
我的 CRS 是 +proj=longlat +datum=WGS84 +no_defs，根据 documentation of raster::buffer()，它不一定将米作为其默认地图单位。 st_crs(ncin) 产生以下详细信息：
st_crs(ncin) # does the current crs take meters as length-unit?
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
以上 LENGTHUNIT["metre",1] 是否意味着缓冲区半径将以米为单位正确缩放？继续提取让我...
NGA_spatial <- SpatialPointsDataFrame(
NGA_panel[108:107], proj4string=ncin@crs, NGA_panel) # spatialize dataset (i.e. make SPDF)
NGA_spatial <- NGA_spatial[sample(nrow(NGA_spatial), 100), ] #### LIMIT to 100 random rows FOR TESTING!
# extract values for households, using a circular buffer with radius 20KM:
data <- raster::extract(ncin,
NGA_spatial, # SPDF with centroids for buffer
buffer = 20000, # buffer size, units depend on CRS (meters or degrees?)
fun=mean, # what value to extract
df=TRUE) # create a data frame?
summary(data[sample(ncol(data), 3)]) # summary statistics of 3 random columns
#> X2005.09.16 X2000.02.15 X2001.02.15
#> Min. :28.80 Min. :29.20 Min. :31.20
#> 1st Qu.:29.77 1st Qu.:30.30 1st Qu.:32.17
#> Median :30.70 Median :33.30 Median :34.10
#> Mean :31.21 Mean :32.46 Mean :33.63
#> 3rd Qu.:32.25 3rd Qu.:33.80 3rd Qu.:34.55
#> Max. :37.10 Max. :35.00 Max. :36.00
显然我提取的值存在一些变化，但我不确定这是通过 20000 米的缓冲区（如我所愿）还是 20000 度的缓冲区（这很荒谬）提取的。
有人能告诉我这是否正确地以米为单位提取，是否有办法控制它（例如通过某种方式绘制缓冲区）？非常感谢！
您需要查看WKT的分组：
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
里面有 3 个 *UNIT 块：
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]
和 2 个轴：
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
第一个适用于椭圆体（以米为单位），另外两个是您感兴趣的，因为它们定义了数据的存储方式（和呈现方式）以及以度为单位（与未投影的预期一样） WGS CRS）。因此，您的缓冲区将以度为单位，正如您所说的那样是荒谬的，只要没有任何中断，它将 return 您的所有数据。
要以米为单位工作，您需要将数据重新投影到平面本地 CRS（最好是对您的数据集等距离的 CRS）。
?raster::extract 参数“buffer”的文档指出
If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters).
+proj=longlat +datum=WGS84明明是(?)latitude/longitude，还是用米吧
如下所示，取自 ?extract 中的示例（略有改动）：
library(raster)
r <- raster(ncols=36, nrows=18, vals=1:(18*36))
xy <- cbind(0, seq(20, 80, by=20))
xy
# [,1] [,2]
#[1,] 0 20
#[2,] 0 40
#[3,] 0 60
#[4,] 0 80
extract(r, xy, buffer=1000000)
#[[1]]
#[1] 234 235 270 271
#[[2]]
#[1] 162 163 198 199
#[[3]]
#[1] 89 90 91 92 126 127
#[[4]]
# [1] 13 14 15 16 17 18 19 20 21 22 23 24 51 52 53 54 55 56 57 58
如您所见，当您向极点移动时，您会获得越来越多的值（即缓冲区内的单元格）。这是预期的，因为纵向距离减小，大致与 cos(latitude).