converting lat long

General
,
1

I am trying to project GRS80 coordinates to WGS84 to get proper lat/ long values.
When I use projectRaster function raster layer dimension change and the value also change. how can I get same same number grid point and value that rsater value and dimention doesnt change

dem
class : RasterLayer
dimensions : 4100, 7401, 30344100 (nrow, ncol, ncell)
resolution : 1000, 1000 (x, y)
extent : -3901000, 3500000, -5100000, -1e+06 (xmin, xmax, ymin, ymax)
crs : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
source : apg21r_1_0_0
names : Band_1
values : 0, 26602.31 (min, max)

dem2<- projectRaster(dem, crs = crs("+init=epsg:4326"))

dem2
class : RasterLayer
dimensions : 4736, 8523, 40364928 (nrow, ncol, ncell)
resolution : 0.0103, 0.00891 (x, y)
extent : 85.90868, 173.6956, -46.81898, -4.621222 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : Band_1
values : 0, 21502.74 (min, max)

As you can notice, dem2 has a different cell number, and the value range changed compare to dem

2

Welcome to the forum

I think we need a lot more information.
See FAQ Asking Questions

3

Here is something that is reproducible:

library(terra)

# sample data in terra package
elev_path <- system.file("ex/elev.tif", package="terra")
elev_wgs84 <- rast(elev_path) # in 4326 by default
elev_grs80 <- project(elev_wgs84, "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs")

There is a bit of an answer in the comments here: r - terra::project seems to be changing the values of a raster - Geographic Information Systems Stack Exchange

There is more information in the notes of terra::project().

4

As William mentioned, dimension and value change with projection.
The initial raster got 90 rows and 95 columns; after projection, 126 rows and 103 columns,

Also, value range and mean everything changes, so projection leads to completely wrong values for further calculations.

I want the exact projected lat long for a given value in the initial raster file.

elev_wgs84 <- rast(elev_path)
elev_wgs84
class : SpatRaster
dimensions : 90, 95, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : elev.tif
name : elevation
**min value : 141 **
max value : 547
elev_grs80 <- project(elev_wgs84, "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs")
elev_grs80
class : SpatRaster
dimensions : 126, 103, 1 (nrow, ncol, nlyr)
resolution : 967.4976, 967.4976 (x, y)
extent : -16226457, -16126804, -4784178, -4662273 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
source(s) : memory
name : elevation
**> min value : 141.8710 **
> max value : 544.9891

5

To maintain the same grid dimensions and values when projecting, first reproject the raster with projectRaster, then use resample to align it back to the original grid dimensions and resolution. Ensure the output matches the original extent and resolution for consistency. This approach keeps the number of grid points and values unchanged.

1 Like
6

When you project a raster by specifying the CRS as the second argument, the dimensions of the projected raster are determined "automagically" and generally don't correspond to the original raster. This means that the number of cells is different, and the values are obviously not the same.
(Here I'm using a more suitable CRS for Luxembourg that does not distort the geometry)

library(terra)
# sample data in terra package
elev_path <- system.file("ex/elev.tif", package="terra")
elev_source <- rast(elev_path) # in 4326 by default
# use LUREF / Luxembourg TM projection (EPSG:2169) 
# https://epsg.io/2169
target_crs <- "epsg:2169"
elev_target <- project(elev_source, target_crs)
# compare source and target
print(elev_source)
print(elev_target)
hist(elev_source)
hist(elev_target)
# not same dimensions, somewhat different values, this is to be expecte

The recommended way of using the project function is to provide a target raster template (dimensions/resolution, extent, and CRS). This will keep the same dimensions; however, the values will generally still differ.

# from the terra::project help page
# if (x is a SpatRaster, the preferred approach is for y to be a SpatRaster as well, 
# serving as a template for the geometry (extent and resolution) of the output SpatRaster. 
# Alternatively, you can provide a coordinate reference system (CRS) description.

# define target raster template
target_ext <- project(ext(elev_source), from = crs(elev_source), to = target_crs)
rast_template <- rast(nrows = nrow(elev_source), ncols = ncol(elev_source), extent = target_ext, crs = target_crs)
# project using target raster template
elev_target <- project(elev_source, rast_template)
# compare source and target
print(elev_source)
print(elev_target)
old.par <- par(pty = "s")
plot(values(elev_source), values(elev_target), pch = 16, xlab = "source", ylab = "target")
par(old.par)
# same dimensions, somewhat different values, this is to be expected

To keep the exact same values, you can try going through the vector format (polygons). It works in reprex, and hopefully, it will also work with your real data.

# create target raster by going through polygons
# convert to polygons
poly_source <- as.polygons(elev_source, round = FALSE, aggregate = FALSE, na.rm = FALSE)
# project to target CRS
poly_target <- project(poly_source, target_crs)
# convert projected polygons back to raster
elev_target <- rasterize(poly_target, rast_template, field = names(elev_source)[1])
# compare source and target
print(elev_source)
print(elev_target)
old.par <- par(pty = "s")
plot(values(elev_source), values(elev_target), pch = 16, xlab = "source", ylab = "target")
par(old.par)
print(all.equal(values(elev_source), values(elev_target)))
# same dimensions, same values
1 Like
7

My input data is Australian population grid 2023 in ESRI Grid format
I could resolve the issue by changing the data SPDF and getting back to the grid

"dem is Australian population grid 2023 in ESRI Grid raster data."

PG <- crop(dem, extent(700000, 2100000, -4250000, -3100000))
PG
class : RasterLayer
dimensions : 1150, 1400, 1610000 (nrow, ncol, ncell)
resolution : 1000, 1000 (x, y)
extent : 7e+05, 2100000, -4250000, -3100000 (xmin, xmax, ymin, ymax)
crs : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
source : memory
names : Band_1
values : 0, 26602.31 (min, max)

Then I extracted the coordinates and population data, created a dataframe

df<-coordinates(PG)
pop_v<-values(PG)

ndt<-data.frame(cbind(df, pop_v))

coordinates(ndt) <- c("x", "y")
proj4string(ndt) <- CRS("+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs") # WGS 84
CRS.old <-CRS("+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs")
CRS.new <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 +units=m")

ndt_pg <- spTransform(ndt, CRS.new)
ndt_pg
class : SpatialPointsDataFrame
features : 1610000
extent : 139.2377, 155.6683, -38.66444, -26.96894 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
variables : 1
names : pop_v
min values : 0
max values : 26602.306640625

By doing so, I got SPDF data with desired coordinates and values; values and number of cells stayed the same.
I can use SPDF data for calculation lat long and return to the same raster format or a different one.