geospatial - Unable to set CRS for list of raster maps using functionloop in R

I have a list of raster map bricks created by cropping existing rasters. I also have a csv metadata file which includes the EPSG code for each raster. The crs info within the cropped rasters, inherited from the originals, seems to be working but is incomplete and missing the datum name; datum shows as 'unknown' when exported to tif and imported to GIS. I wish to update the entire list by fetching the individual EPSG codes from the metadata file and using raster:crs or sp:CRS or etc. I can easily do this manually for each raster using a variety of commands, but the operation fails consistently when I try to loop it over the entire list.

This is what I want to achieve:

ORIGINAL:

crs(cropped_rasterdata[[1]])
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on GRS 1980 ellipsoid",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1],
ID["EPSG",7019]]],
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]]]]

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

NEW:

crs(cropped_rasterdata[[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]]

However, when I try this:

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

I get no result:

crs(cropped_rasterdata[[2]])
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on GRS 1980 ellipsoid",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1],
ID["EPSG",7019]]],
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]]]]

No errors are shown.

The aim is to get a function working and then loop it down the list with something like:

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

Needs to be applicable to raster lists of varying length.

Any assistance with this would be much appreciated.

EDIT: Reproducible Example

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
# ?

This is what I want to achieve:

ORIGINAL:

crs(cropped_rasterdata[[1]])
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on GRS 1980 ellipsoid",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1],
ID["EPSG",7019]]],
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]]]]

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

NEW:

crs(cropped_rasterdata[[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]]

However, when I try this:

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

I get no result:

crs(cropped_rasterdata[[2]])
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=GRS80 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on GRS 1980 ellipsoid",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1],
ID["EPSG",7019]]],
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]]]]

No errors are shown.

The aim is to get a function working and then loop it down the list with something like:

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

Needs to be applicable to raster lists of varying length.

Any assistance with this would be much appreciated.

EDIT: Reproducible Example

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
# ?

Share Improve this question edited Nov 19, 2024 at 12:29 asked Nov 19, 2024 at 5:44 rommel 133 bronze badges

At least, there is a typo in your fun_crs; a ) is missing. Btw, without the option to specify different objects, this function appears useless to me. – Friede Commented Nov 19, 2024 at 9:58
If cropped_rasterdata and raster_metadata$EPSG are of same length you can probably do sth along the lines: cropped_rasterdata = lapply(seq_along(cropped_rasterdata), \(x) crs(cropped_rasterdata[[x]]) = raster_metadata$EPSG[x]). – Friede Commented Nov 19, 2024 at 10:02
Please use spaces, I edited again. It makes code easier to read. Also you night to assign for attempts as I showed in my previous comment. Side note: Do not use for (i in 1:nrow(>?<)), use for (i in seq(nrow(>?<))). – Friede Commented Nov 19, 2024 at 10:24
@Friede Thanks, I edited the typo. I’m not sure what you mean by useless? I’ve shown how the same command works fine with one item in a list. Admittedly I’m no R wizard but I’m perplexed as to why, when I replace the constants (in this case row numbers) with variables in a function, the very same command does not do the same thing. The raster list and the metadata are the same length and the information is in the same row. Thanks for your suggestion, I’ll give it a try. – rommel Commented Nov 19, 2024 at 10:30
1 Yes, the EPSGs are the same in this case, but obviously not necessarily. Your answer does help me; works exactly as I wanted it to, cheers. I still don't understand why my useless function doesn't work though when the underlying command works perfectly. – rommel Commented Nov 19, 2024 at 13:39

| Show 5 more comments

1 Answer 1

Sorted by: Reset to default 0

Before

> 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]]]]

I suspect

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

works. After:

> 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]]

科技改变生活-雨落星辰 - 所有的伟大,都源于一个勇敢的开始

geospatial - Unable to set CRS for list of raster maps using functionloop in R - Stack Overflow

1 Answer 1

与本文相关的文章

评论列表(0)