r - How to change coordinate system of extent to match raster? - Stack Overflow

admin管理员组

文章数量:1435859

Here is my current R code:

library(ebirdst)
library(raster)
library(sf)
library(terra)

# specify parameters
species = "American Robin"
region_extent <- ext(-76.353957, -75.246633, 44.965633, 45.536983)

# get occurrence data
ebirdst_download_status(species, download_abundance = FALSE, download_occurrence = TRUE, force = FALSE)
occurrence_raster <- load_raster(species, product=c("occurrence"))

# filter to county
occurrence_raster <- project(occurrence_raster, crs(region_extent))
cropped_raster <- crop(occurrence_raster, region_extent)

This works, but projecting the raster onto the coordinate system of the extent in the second last step is very slow. Is there a way I can change the coordinate system of the extent instead? When I print the occurrence_raster, here's what I get:

class       : SpatRaster 
dimensions  : 5630, 13511, 52  (nrow, ncol, nlyr)
resolution  : 2962.807, 2962.809  (x, y)
extent      : -20015109, 20015371, -6673060, 10007555  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 

Here is my current R code:

library(ebirdst)
library(raster)
library(sf)
library(terra)

# specify parameters
species = "American Robin"
region_extent <- ext(-76.353957, -75.246633, 44.965633, 45.536983)

# get occurrence data
ebirdst_download_status(species, download_abundance = FALSE, download_occurrence = TRUE, force = FALSE)
occurrence_raster <- load_raster(species, product=c("occurrence"))

# filter to county
occurrence_raster <- project(occurrence_raster, crs(region_extent))
cropped_raster <- crop(occurrence_raster, region_extent)

This works, but projecting the raster onto the coordinate system of the extent in the second last step is very slow. Is there a way I can change the coordinate system of the extent instead? When I print the occurrence_raster, here's what I get:

class       : SpatRaster 
dimensions  : 5630, 13511, 52  (nrow, ncol, nlyr)
resolution  : 2962.807, 2962.809  (x, y)
extent      : -20015109, 20015371, -6673060, 10007555  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 

• r
• raster

Share Improve this question Follow asked Nov 15, 2024 at 19:32 Jan HuusJan Huus 1031010 bronze badges

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 2

You can:

1. Create a cropping object with a defined CRS (assumed here to be WGS84/EPSG:4326)
2. Get CRS of the SpatRaster, and use to project the cropping object
3. Use projected cropping object to crop() the SpatRaster

Note that crop() works by returning rows/columns so the extent of the projected cropping object and the cropped SpatRaster will differ in this example.

library(ebirdst)
library(raster)
library(sf)
library(terra)

# Specify parameters
species = "American Robin"

# Set crop extent as SpatVector with defined CRS
region_extent <- vect(ext(-76.353957, -75.246633, 44.965633, 45.536983), 
                      crs = "EPSG:4326")

# Get occurrence data
ebirdst_download_status(species,
                        download_abundance = FALSE,
                        download_occurrence = TRUE,
                        force = FALSE)

occurrence_raster <- load_raster(species, product = c("occurrence"))

# Get CRS from occurrence_raster
prj_crs <- crs(occurrence_raster)

# Project region to match CRS of occurrence_raster
region_prj <- project(region_extent, prj_crs)

region_prj
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 0  (geometries, attributes)
# extent      : -6007065, -5860692, 4999956, 5063487  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 

# Crop occurrence_raster to projected region extent
cropped_raster <- crop(occurrence_raster, region_prj)

cropped_raster
# class       : SpatRaster 
# dimensions  : 21, 49, 52  (nrow, ncol, nlyr)
# resolution  : 2962.807, 2962.809  (x, y)
# extent      : -6006959, -5861782, 5000408, 5062626  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
# source(s)   : memory
# varname     : amerob_occurrence_median_3km_2022 
# names       : 2022-01-04, 2022-01-11, 2022-01-18, 2022-01-25, 2022-02-01, 2022-02-08, ... 
# min values  :  0.0000000,  0.0000000,  0.0000000,  0.0000000,   0.000000,  0.0000000, ... 
# max values  :  0.7491944,  0.7453837,  0.7386807,  0.7669981,   0.759219,  0.7370626, ...

本文标签： rHow to change coordinate system of extent to match rasterStack Overflow