无法使用 R 中的函数/循环为栅格地图列表设置 CRS

Question

我有一个通过裁剪现有栅格创建的栅格地图砖块列表。 我还有一个 csv 元数据文件，其中包含每个栅格的 EPSG 代码。 裁剪后的栅格中的 crs 信息继承自原始栅格，似乎有效，但不完整且缺少基准名称；导出到 tif 并导入到 GIS 时，数据显示为“未知”。 我希望通过从元数据文件中获取各个 EPSG 代码并使用 raster:crs 或 sp:CRS 等来更新整个列表。我可以使用各种命令轻松地为每个栅格手动执行此操作，但操作始终会失败我尝试将其循环到整个列表。

这就是我想要实现的目标：

原文：

crs(cropped_rasterdata[[1]])
坐标参考系：
已弃用的 Proj.4 表示形式：+proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 代表：
GEOGCRS[“未知”，
DATUM[“基于 GRS 1980 椭球的未知”，
椭圆体["GRS 1980",6378137,298.257222101,
LENGTHUNIT["米",1],
ID["EPSG",7019]]],
PRIMEM["格林威治",0,
ANGLEUNIT["度",0.0174532925199433],
ID["EPSG",8901]],
CS[椭球,2],
AXIS["经度",东,
订单[1]，
ANGLEUNIT["度",0.0174532925199433,
ID["EPSG",9122]]],
轴[“纬度”，北，
订单[2]，
ANGLEUNIT["度",0.0174532925199433,
ID["EPSG",9122]]]]

crs(cropped_rasterdata[[1]]) <- raster_metadata$EPSG[1]

新：

crs(cropped_rasterdata[[1]])
坐标参考系：
已弃用的 Proj.4 表示形式：+proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 代表：
GEOGCRS["GDA2020",
DATUM[“2020 年澳大利亚地心基准面”，
椭圆体["GRS 1980",6378137,298.257222101,
LENGTHUNIT["米",1]]],
PRIMEM["格林威治",0,
ANGLEUNIT["度",0.0174532925199433]],
CS[椭球,2],
AXIS["大地纬度 (Lat)",北,
订单[1]，
ANGLEUNIT["度",0.0174532925199433]],
AXIS["大地经度 (Lon)",东,
订单[2]，
ANGLEUNIT["度",0.0174532925199433]],
使用方法[
范围[“大地测量、地籍、工程测量、地形测绘。”],
区域[“澳大利亚包括豪勋爵岛、麦觉理岛、阿什莫尔和卡地亚群岛、圣诞岛、科科斯（基林）群岛、诺福克岛。所有陆上和近海。”],
BBOX[-60.55,93.41,-8.47,173.34]],
ID["EPSG",7844]]

但是，当我尝试这个时：

fun_crs <- function(x) { crs(cropped_rasterdata[[x]]) <- raster_metadata$EPSG[x] }
fun_crs(2)

我没有得到结果：

crs(cropped_rasterdata[[2]])
坐标参考系：
已弃用的 Proj.4 表示形式：+proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 代表：
GEOGCRS[“未知”，
DATUM[“基于 GRS 1980 椭球的未知”，
椭圆体["GRS 1980",6378137,298.257222101,
LENGTHUNIT["米",1],
ID["EPSG",7019]]],
PRIMEM["格林威治",0,
ANGLEUNIT["度",0.0174532925199433],
ID[“EPSG”，8901]]，
CS[椭球,2],
AXIS["经度",东,
订单[1]，
ANGLEUNIT["度",0.0174532925199433,
ID["EPSG",9122]]],
轴[“纬度”，北，
订单[2]，
ANGLEUNIT["度",0.0174532925199433,
ID["EPSG",9122]]]]

没有显示错误。

目标是让一个函数正常工作，然后将其在列表中循环，如下所示：

for(i in 1:nrow) fun_crs(i)

需要适用于不同长度的栅格列表。

如有任何帮助，我们将不胜感激。

编辑：可重现的示例

library (raster)
library (sp)
  
# create raster bricks
brick1 <- brick()
brick2 <- brick()
brick3 <- brick()
  
# create brick list
brick_list <- list(brick1, brick2, brick3)

# create EPSG list
EPSG <- rep(7844,3)
DF <- data.frame(EPSG)

# check brick1 crs
crs(brick_list[[1]])

# update brick1 crs
crs(brick_list[[1]]) <- DF$EPSG[1]

# check brick1 crs (this works)
crs(brick_list[[1]])

# function to apply in loop
fun_crs <- function(x) { crs(brick_list[[x]]) <- DF$EPSG[x] }

# test function on brick2
fun_crs(2)

# check brick2 crs (this did not work)
crs(brick_list[[2]])

# apply/loop function to update all bricks in list
# ?

Answer 1

之前

> lapply(brick_list, crs)
[[1]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
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]]]] 

[[2]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
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]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
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]]]] 

我怀疑

brick_list = lapply(seq_along(brick_list), \(i) {
  crs(brick_list[[i]]) = DF$EPSG[i]
  brick_list[[i]] })

有效。之后：

> lapply(brick_list, crs)
[[1]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs 
WKT2 2019 representation:
GEOGCRS["GDA2020",
    DATUM["Geocentric Datum of Australia 2020",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Geodesy, cadastre, engineering survey, topographic mapping."],
        AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
        BBOX[-60.55,93.41,-8.47,173.34]],
    ID["EPSG",7844]] 

[[2]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs 
WKT2 2019 representation:
GEOGCRS["GDA2020",
    DATUM["Geocentric Datum of Australia 2020",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Geodesy, cadastre, engineering survey, topographic mapping."],
        AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
        BBOX[-60.55,93.41,-8.47,173.34]],
    ID["EPSG",7844]] 

[[3]]
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs 
WKT2 2019 representation:
GEOGCRS["GDA2020",
    DATUM["Geocentric Datum of Australia 2020",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Geodesy, cadastre, engineering survey, topographic mapping."],
        AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
        BBOX[-60.55,93.41,-8.47,173.34]],
    ID["EPSG",7844]]