Reputation: 1215
Converting the Coordinate Reference System (CRS) of a Stack of 70+ Raster Files in stack::raster () Format
Issue
I have a batch of 70+ raster files with the object name 'ncin_SST' that I created using the function stack::raster() in the raster package. I want to change the coordinate reference system (CRS) WG84/34N.
My ultimate aim is to interpolate the point data (dolphin IDs) in a shapefile to obtain the sea surface temperature (SST) values from each individual raster file within the stack of 70+ rasters.
The source of the raster files are multiple Aqua Modis netCDF files contained in one folder in my documents and I downloaded the data from the 'Ocean Color Project' from NASA. My goal is to extract the average SST per ID over the time period across all 70+ raster files from 2016-to 2021.
Firstly, I need to ensure that the shapefile and raster files have the same CRS and I have tried various different methods to achieve this output; however, unsuccessfully.
Unfortunately, I cannot share my data online due to ownership of the data.
Many thanks in advance if anyone can help with this error message.
R-code:
##Open all 70+ Aqua Modis netCDF files from their folder
filenames = list.files('~/Ocean_ColorSST_2016_2021', pattern='*.nc',full.names=TRUE)
##Stack all the netCDF files and extract the variable "sst"
ncin_SST <- raster::stack(filenames, varname = "sst")
Print Results
print(ncin_SST)
class : RasterStack
dimensions : 4320, 8640, 37324800, 59 (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -180, 180, -90.00001, 90 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
names : Sea.Surface.Temperature.1, Sea.Surface.Temperature.2, Sea.Surface.Temperature.3, Sea.Surface.Temperature.4, Sea.Surface.Temperature.5, Sea.Surface.Temperature.6, Sea.Surface.Temperature.7, Sea.Surface.Temperature.8, Sea.Surface.Temperature.9, Sea.Surface.Temperature.10, Sea.Surface.Temperature.11, Sea.Surface.Temperature.12, Sea.Surface.Temperature.13, Sea.Surface.Temperature.14, Sea.Surface.Temperature.15, ...
Here's a printout of one single AQUA MODIS FILE:
File AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc (NC_FORMAT_NETCDF4):
3 variables (excluding dimension variables):
short sst[lon,lat] (Chunking: [87,44]) (Compression: shuffle,level 4)
long_name: Sea Surface Temperature
scale_factor: 0.00499999988824129
add_offset: 0
units: degree_C
standard_name: sea_surface_temperature
_FillValue: -32767
valid_min: -1000
valid_max: 10000
display_scale: linear
display_min: -2
display_max: 45
unsigned byte qual_sst[lon,lat] (Chunking: [87,44]) (Compression: shuffle,level 4)
long_name: Quality Levels, Sea Surface Temperature
_FillValue: 255
valid_min: 0
valid_max: 5
unsigned byte palette[eightbitcolor,rgb] (Contiguous storage)
4 dimensions:
lat Size:4320
long_name: Latitude
units: degrees_north
standard_name: latitude
_FillValue: -999
valid_min: -90
valid_max: 90
lon Size:8640
long_name: Longitude
units: degrees_east
standard_name: longitude
_FillValue: -999
valid_min: -180
valid_max: 180
rgb Size:3 (no dimvar)
eightbitcolor Size:256 (no dimvar)
[1] ">>>> WARNING <<< attribute data_bins is an 8-byte value, but R"
[1] "does not support this data type. I am returning a double precision"
[1] "floating point, but you must be aware that this could lose precision!"
61 global attributes:
product_name: AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc
instrument: MODIS
title: MODISA Level-3 Standard Mapped Image
project: Ocean Biology Processing Group (NASA/GSFC/OBPG)
platform: Aqua
temporal_range: month
processing_version: R2019.0
date_created: 2019-12-17T18:25:34.000Z
history: l3mapgen par=AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc.param
l2_flag_names: LAND,HISOLZEN
time_coverage_start: 2016-09-01T00:00:00.000Z
time_coverage_end: 2016-10-01T02:59:59.000Z
start_orbit_number: 76217
end_orbit_number: 76655
map_projection: Equidistant Cylindrical
latitude_units: degrees_north
longitude_units: degrees_east
northernmost_latitude: 90
southernmost_latitude: -90
westernmost_longitude: -180
easternmost_longitude: 180
geospatial_lat_max: 90
geospatial_lat_min: -90
geospatial_lon_max: 180
geospatial_lon_min: -180
latitude_step: 0.0416666679084301
longitude_step: 0.0416666679084301
sw_point_latitude: -89.9791641235352
sw_point_longitude: -179.97917175293
spatialResolution: 4.64 km
geospatial_lon_resolution: 0.0416666679084301
geospatial_lat_resolution: 0.0416666679084301
geospatial_lat_units: degrees_north
geospatial_lon_units: degrees_east
number_of_lines: 4320
number_of_columns: 8640
measure: Mean
suggested_image_scaling_minimum: -2
suggested_image_scaling_maximum: 45
suggested_image_scaling_type: LINEAR
suggested_image_scaling_applied: No
_lastModified: 2019-12-17T18:25:34.000Z
Conventions: CF-1.6 ACDD-1.3
institution: NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group
standard_name_vocabulary: CF Standard Name Table v36
naming_authority: gov.nasa.gsfc.sci.oceandata
id: AQUA_MODIS.20160901_20160930.L3b.MO.SST.nc/L3/AQUA_MODIS.20160901_20160930.L3b.MO.SST.nc
license: https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
creator_name: NASA/GSFC/OBPG
publisher_name: NASA/GSFC/OBPG
creator_email: [email protected]
publisher_email: [email protected]
creator_url: https://oceandata.sci.gsfc.nasa.gov
publisher_url: https://oceandata.sci.gsfc.nasa.gov
processing_level: L3 Mapped
cdm_data_type: grid
keywords: Earth Science > Oceans > Ocean Optics > Sea Surface Temperature
keywords_vocabulary: NASA Global Change Master Directory (GCMD) Science Keywords
data_bins: 20740391
data_minimum: -1.80000019073486
data_maximum: 39.9599990844727
Change the CRS to WG84/UTM 34
ncin_SST_34 <- projectRaster(ncin_SST, crs = "+proj=utm + zone=34 + datum=WGS84")
Error Message
Error in rgdal::checkCRSArgs_ng(uprojargs = uprojargs, SRS_string = SRS_string, :
length(SRS_string) == 1L is not TRUE
In addition: Warning messages:
1: In if (is.na(x)) { :
the condition has length > 1 and only the first element will be used
2: In if (x == "") { :
the condition has length > 1 and only the first element will be used
3: In if (substr(x, 1, 1) == "+") { :
the condition has length > 1 and only the first element will be used
Upvotes: 0
Views: 634
Answers (1)
Reputation: 47051
To open a ncdf file with raster data you can simply do
library(terra)
x <- rast("filename.nc")
if there are multiple variables per file you can select one like this
x <- rast("filename.nc", "SST")
(or open all of them as a SpatRasterDataset with sds("filename.nc")
If you have a vector ff that contains multiple filenames of ncdf files for the same area (same extent and resolution) you can do
x <- rast(ff, "SST")
While it is possible to transform your data to another coordinate reference system such as UTM like this
p <- project(x, "+proj=utm + zone=34 + datum=WGS84")
That is in most cases not a good idea. First, you should really provide a template raster to project to (see ?project), but more importantly, projecting raster data distorts the values. So you should avoid that if possible in data analysis (it is ok to make a map in most cases). Instead project the spatial vector data to the CRS of the raster data.
Upvotes: 2
Related Questions
- Writing a raster stack to file with changed names
- How can I extract data from a raster stack based on a list of lat long?
- Select specific rasters from stack and make a new stack in R
- Identify CRS in raster file
- Conversion of Coordinate reference system in Raster objects in R
- Converting raster stack extent from meters to lat/lon coordinates
- Naming layers in a raster stack and extracting data
- Change coordinates to a series of rasters with multiple layers in R
- how to combine raster layers from a stack into a data frame R
- stack raster images from a folder in R