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