Reprojected raster data shifts East
I'm working in R using terra and sf. I am working with some data products that are in sinusoidal projection (see crs string below, I can't change the starting product projection) and I need to reproject them into a different equal area projection (currently selected EPSG:9311, but open to suggestions) so that other people can use and interpret the maps more easily. The issue is, when I reproject, the data seem to shift eastward about 2 grid cells. I think this has to do with some kind of rounding bias, but I really don't have an explanation of why this is occurring. Can you help?
library(terra)
library(sf)
library(dplyr)
library(rnaturalearth)
# get sinusoidal projection and resolution of original raster
a_crs <- 'PROJCRS["unnamed",
BASEGEOGCRS["unnamed ellipse",
DATUM["unknown",
ELLIPSOID["unnamed", 6371007.181, 0,
LENGTHUNIT["metre", 1,
ID["EPSG", 9001]]]],
PRIMEM["Greenwich", 0,
ANGLEUNIT["degree", 0.0174532925199433,
ID["EPSG", 9122]]]],
CONVERSION["Sinusoidal",
METHOD["Sinusoidal"],
PARAMETER["Longitude of natural origin", 0,
ANGLEUNIT["degree", 0.0174532925199433],
ID["EPSG", 8802]],
PARAMETER["False easting", 0,
LENGTHUNIT["metre", 1],
ID["EPSG", 8806]],
PARAMETER["False northing", 0,
LENGTHUNIT["metre", 1],
ID["EPSG", 8807]]],
CS[Cartesian, 2],
AXIS["easting", east,
ORDER[1],
LENGTHUNIT["metre", 1,
ID["EPSG", 9001]]],
AXIS["northing", north,
ORDER[2],
LENGTHUNIT["metre", 1,
ID["EPSG", 9001]]]]'
a_res <- c(2962.807, 2962.8090)
#Get the state of Connecticut (the spatial error is more obvious in smaller states)
state <- ne_states(country = "United States of America", returnclass = "sf") |>
filter(name == 'Connecticut') |>
st_transform(a_crs)
# Define the extent of the test grid (xmin, xmax, ymin, ymax)
extent_grid <- ext(state)
# Create the raster with the specified extent and resolution
raster_grid <- rast(extent = extent_grid, resolution = a_res) # ~3km resolution
# Set the CRS of the raster
crs(raster_grid) <- a_crs
# assign a grid of values to the test grid
values(raster_grid) <- c(1,rep(0,9))
# grid cell touches the state or has centroid within the state
mask <- rasterize(vect(state), raster_grid, background=NA, touches = T)
raster_grid[!is.na(mask)] <- raster_grid[!is.na(mask)] + 3
# grid cell has centroid within the state
mask <- rasterize(vect(state), raster_grid, background=NA, fun = 'sum')
raster_grid[!is.na(mask)] <- raster_grid[!is.na(mask)] + 3
# view the starting (terrible) data projection and alignment with state vector geometry
plot(raster_grid)
plot(state$geometry, add = T)
# transform the state and raster to the new crs
transformed_state <- st_transform(state, crs('EPSG:9311'))
transformed_raster <- project(raster_grid, crs('EPSG:9311'), res = a_res, method='near')
# view problematic reprojection
plot(transformed_raster) #Why??!
plot(transformed_state$geometry, add = T)
This is the expected/desired output
# Get the current extent
current_extent <- ext(transformed_raster)
# Calculate the new extent by shifting 2 cells left
new_extent <- c(
xmin = current_extent$xmin - 2 * a_res[1],
xmax = current_extent$xmax - 2 * a_res[1],
ymin = current_extent$ymin,
ymax = current_extent$ymax)
# Create a new raster with the adjusted extent
shifted_raster <- rast(ncol= ncol(transformed_raster), nrow= nrow(transformed_raster),
ext= new_extent, crs= crs(transformed_raster))
# Copy the values from the original raster to the new raster
values(shifted_raster) <- values(transformed_raster)
# Plot shifted alignment
plot(shifted_raster)
plot(transformed_state$geometry, add = T)
This also occurs with other reprojection methods, such as bilinear interpolation. I wasn't sure if there was a trigonometric reason for the shift, like perhaps it was sqrt(2)/2 in the south-easterly direction, but that didn't seem to improve the accuracy either.
EDIT: I also ran this in the terminal using gdalwarp and still get the same issue when projecting to EPSG:9311
gdalwarp -wo SAMPLE_GRID=YES -t_srs EPSG:9311 -r near raster_grid.tif proj_raster_grid_CT.tif
Your issue is that you are selecting the wrong CRS to project to. In your case, your data are outside the projected bounds of NAD27 US National Atlas Equal Area/EPSG:9311. Data outside the projected bounds of a CRS will be subjected to increased distortions like those you are experiencing. The further they are from the projected bounds, the greater the distortion.
First, an illustration of your data compared the the projected bounds of EPSG:9311:
library(terra)
library(sf)
library(dplyr)
library(rnaturalearth)
library(ggplot2)
# Project ne_countries()
world9311 <- st_transform(ne_countries(), 9311) |>
select()
# EPSG:9311 project bounds
ext9311 <- st_as_sf(vect(ext(-141.01, -44.0, 40.0, 83.17), crs = "EPSG:4326")) |>
st_transform(9311)
# Connecticut EPSG:9311 sf
state9311 <- ne_states(country = "United States of America", returnclass = "sf") |>
filter(name == 'Connecticut') |>
st_transform(9311)
# Note that the latitudinal extent of state9311 falls outside EPSG:9311
ext(ext9311)
# SpatExtent : -3317821.01925665, 4316225.9634226, 302103.872897179, 4462577.97807183
# (xmin, xmax, ymin, ymax)
ext(state9311)
# SpatExtent : 2142663.86101191, 2295220.74938045, -92217.4826628708, 67242.1773451936
# (xmin, xmax, ymin, ymax)
ggplot() +
geom_sf(data = world9311, colour = "grey75") +
geom_sf(data = ext9311, aes(colour = "EPSG:9311\nProjected Bounds"),
fill = NA) +
geom_sf(data = state9311, aes(fill = "Connecticut")) +
scale_colour_manual(name = "",
values = "#DF536B") +
scale_fill_manual(name = "",
values = "#490092") +
guides(color = guide_legend(override.aes = list(fill = "white"))) +
coord_sf(xlim = c(ext(ext9311)[1], ext(ext9311)[2]),
ylim = c(ext(state9311)[3], ext(ext9311)[4])) +
theme(legend.position = "bottom")
An example where your data are within the projected bounds of EPSG:9311 using Qikiqtaq as a proxy for Connecticut. As per my suggestion in a previous edit of this answer, this example uses a template for project(), which is the preferred method for projecting a SpatRaster using terra:
# Example data
Qikiqtaq <- st_cast(world9311, "POLYGON") |>
slice(33) |>
st_transform(a_crs)
# Create raster with the specified extent and resolution
raster_grid <- rast(ext(Qikiqtaq), resolution = a_res, crs = a_crs)
# Assign a grid of values to raster_grid
values(raster_grid) <- c(1, rep(0, 9))
# grid cell touches the state or has centroid within the state
mask <- rasterize(vect(Qikiqtaq), raster_grid, background=NA, touches = T)
raster_grid[!is.na(mask)] <- raster_grid[!is.na(mask)] + 3
# grid cell has centroid within the state
mask <- rasterize(vect(Qikiqtaq), raster_grid, background=NA, fun = 'sum')
raster_grid[!is.na(mask)] <- raster_grid[!is.na(mask)] + 3
# Project state to EPSG:9311
transformed_state <- st_transform(Qikiqtaq, 9311)
# Create SpatRaster to use as project() template
ext9311 <- rast(ext(transformed_state),
res = a_res,
crs = "EPSG:9311")
# Project raster_grid to EPSG:9311
transformed_raster <- project(raster_grid,
ext9311,
method = "near")
plot(transformed_raster)
plot(transformed_state$geometry, add = TRUE)
Note that the vector and raster are much more coincident:
Finally, to correct your issue, you will need to choose a CRS that encompasses your data. This example uses North America Albers Equal Area Conic/ESRI:102008, but double-check if there is a more appropriate CRS for your actual data:
# Project state to North America Albers Equal Area Conic/ESRI:102008
transformed_state <- st_transform(state, "ESRI:102008")
# Create SpatRaster to use as project() template
ext102008 <- rast(ext(transformed_state),
res = a_res,
crs = "ESRI:102008")
# Project raster_grid to North America Albers Equal Area Conic/ESRI:102008
transformed_raster <- project(raster_grid,
ext102008,
method = "near")
plot(transformed_raster)
plot(transformed_state$geometry, add = TRUE)
- rlanguageserver: Disable specific linter and avoid line-breaking
- Replace truncated word in a string
- Add label to right y axis in R ggplot
- Popup offset without using markers
- Multi marker depending a variable in a data frame
- Drawing numerous polygons
- Select icons or CircleMarkers based on factor variable
- Why is a map not displayed when I call leaflet in a function
- Display SpatialPolygonsDataFrame on leaflet map
- Leaflet in R not showing variable circle colors
- Overlay two population pyramid plots in one ggplot graph
- How to most efficiently filter a dataframe conditionally of values in another one, in the tidyverse framework?
- R group_by multiple conditionnal aggregation functions dplyr
- Conditionally split strings cells in a dataframe R
- How to have Shiny app render multiple plots
- Disable Shiny button conditionally to a rHandsontable format
- Is there a way to render dynamic headings in Quarto using R?
- What is the time complexity of sample?
- Add geom_text with x bar and round values
- Color-fill ggplot2 maps with a single color
- lapply a ggplot function on multiple inputs with consistent groups
- Is there a way to save text in R without it actively implementing \n?
- How to get number of unique observations using tbl_summary
- Using if_else in a mutate command yields an error about false having the wrong size
- Rotate x-axis value in dendogram using fviz_dend
- Export data frames in list to xlsx with named sheets
- Query all primary roads from PostGis Database filled with OSM Data within a boundary
- Write a list in r to an excel file
- Calculating net assets with missing values: discrepancy in results
- Error in curl::curl_fetch_memory(url, handle = handle) : Couldn't connect to server in R: oauth2.0_token())