I have raster data in the geographic coordinate system EPSG:4326 with units in degrees. I would like to convert my raster to a projected coordinate system with units in meters. I chose +proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs, but I am open to other projection systems as long as the units are in meters. I am getting unexpected results. I expected the cells to have the same values as in EPSG:4326.

Here a figure:

Here is the dataset: .csv?rlkey=p6raz1q0gxjujyt9cp8ilqw33&st=a0b2rufg&dl=0

Here is the code:

data_df <- read.csv("C:/Users/Sophie/Downloads/Test.csv")
data_r <- tidyterra::as_spatraster(tibble::as_tibble(data_df, xy = TRUE), xycols = 1:2, crs = "EPSG:4326")
plot(data_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

test_r <- terra::project(data_r, "+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs")
plot(test_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

I have raster data in the geographic coordinate system EPSG:4326 with units in degrees. I would like to convert my raster to a projected coordinate system with units in meters. I chose +proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs, but I am open to other projection systems as long as the units are in meters. I am getting unexpected results. I expected the cells to have the same values as in EPSG:4326.

Here a figure:

Here is the dataset: https://www.dropbox/scl/fi/i1d5k12x89nrkmihdf6t4/Test.csv?rlkey=p6raz1q0gxjujyt9cp8ilqw33&st=a0b2rufg&dl=0

Here is the code:

data_df <- read.csv("C:/Users/Sophie/Downloads/Test.csv")
data_r <- tidyterra::as_spatraster(tibble::as_tibble(data_df, xy = TRUE), xycols = 1:2, crs = "EPSG:4326")
plot(data_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

test_r <- terra::project(data_r, "+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs")
plot(test_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

• r
• raster
• projection
• terra

Share Improve this question Follow asked Mar 8 at 0:54 Sophie PèreSophie Père 6111 silver badge66 bronze badges

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 1

What you see happens because you are using the default "method="bilinear".

If the raster cell values represent categories that should not be averaged, you may want to use "method="near" instead. There are other methods available, see ?terra::project.

Your example data

library(terra)
d <- read.csv("../Downloads/Test.csv")
r <- rast(d, type="xyz", crs="EPSG:4326")
to <- "+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m"

Solution

p <- project(r, to, method="near")
plot(p)

But since the cell values are fractions the result you were having may actually be better. You could also use "method=average".

---

I wonder why you want the units to be meter. You are mistaken if you think this is useful to compute areas.