最新消息:Welcome to the puzzle paradise for programmers! Here, a well-designed puzzle awaits you. From code logic puzzles to algorithmic challenges, each level is closely centered on the programmer's expertise and skills. Whether you're a novice programmer or an experienced tech guru, you'll find your own challenges on this site. In the process of solving puzzles, you can not only exercise your thinking skills, but also deepen your understanding and application of programming knowledge. Come to start this puzzle journey full of wisdom and challenges, with many programmers to compete with each other and show your programming wisdom! Translated with DeepL.com (free version)

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

matteradmin7PV0评论

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 
Share Improve this question asked Nov 15, 2024 at 19:32 Jan HuusJan Huus 10310 bronze badges
Add a comment  | 

1 Answer 1

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, ...
Post a comment

comment list (0)

  1. No comments so far