Using sp, terra, raster packages in R to convert Albers Equal Area to Lat/Long

Question

I have a .tif file about the ground temperature of a certain region on Earth and I would like to find out the coordinates of the region.

Here is a link to the file I am working on: 2017001D.tif

Using raster package in R, I was able to load the .tif file as RasterLayer class.

I can extract the coordinate information as

> T001D = raster::raster("2017001D.tif")
> T001D
class      : RasterLayer 
dimensions : 4255, 5213, 22181315  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : -2926932, 2286068, 1740497, 5995497  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : 2017001D.tif 
names      : X2017001D 

> coords <- raster::xyFromCell(T001D, seq_len(ncell(T001D)))
> head(coords)
            x       y
[1,] -2926432 5994997
[2,] -2925432 5994997
[3,] -2924432 5994997
[4,] -2923432 5994997
[5,] -2922432 5994997
[6,] -2921432 5994997

I also used terra package to load it as SpatRaster class and when I do so I found more info:

> T001D = terra::rast("2017001D.tif")
> T001D
class       : SpatRaster 
dimensions  : 4255, 5213, 1  (nrow, ncol, nlyr)
resolution  : 1000, 1000  (x, y)
extent      : -2926932, 2286068, 1740497, 5995497  (xmin, xmax, ymin, ymax)
coord. ref. : AEA_WGS_1984 
source      : 2017001D.tif 
name        : 2017001D 

The coordinate system info seems to be AEA_WGS_1984. I looked for it online and found this post:How to convert WGS84 to Lat/Long using R which is similar to my question except I don't have a "zone" number.

It mentioned sp package and I feel like I need help with the functions in it now, such as the CRS syntax in spTransform() function. Can someone help me with this? Thank you

Answer 1

You can do

library(terra)
T001D <- rast("2017001D.tif")
coords <- crds(T001D, na.rm=FALSE)
geo_coords <- project(coords, crs(T001D), "+proj=longlat")

But why would you do that? You probably are approaching what you are after the wrong way. A more common approach (but necessarily what you should do either would be to use project

pT001D <- project(T001D, "+proj=longlat")

That can change the number of rows a bit and columns because it is a different raster geometry; although you can control that. See ?project

Answer 2

You can use terra package to have latitude and longitude in degree decimals like

library(terra)

T001D = terra::rast("2017001D.tif")

T001D_latlon = project(x = T001D, y = "epsg:4326")

coords <- xyFromCell(T001D_latlon, seq_len(ncell(T001D_latlon)))
head(coords)