将 70 多个光栅文件的坐标参考系统 (CRS) 转换为 stack::raster () 格式

Converting the Coordinate Reference System (CRS) of a Stack of 70+ Raster Files in stack::raster () Format

问题

我有一批 70+ 光栅文件,对象名称为 'ncin_SST',我使用函数 stack::raster()raster 包 中。我想更改 坐标参考系统 (CRS) WG84/34N

我的最终目标是在 shapefile 中插入点数据(海豚 ID)以获取 海面温度 (SST) 来自堆栈中每个单独栅格文件的值70 多个栅格。

光栅文件的来源是 多个 Aqua Modis netCDF 文件 包含在我文档的一个文件夹中,我从 'Ocean Color Project' 来自 NASA。我的目标是从 2016 年到 2021 年的所有 70 多个光栅文件中提取每个 ID 的平均 SST。

首先,我需要确保 shapefile 和栅格文件具有相同的 CRS,并且我尝试了各种不同的方法来实现此输出;然而,没有成功。

很遗憾,由于数据的所有权,我无法在线共享我的数据。

非常感谢任何人可以帮助解决此错误消息。

R-代码:

##Open all 70+ Aqua Modis netCDF files from their folder

filenames = list.files('~/Ocean_ColorSST_2016_2021', pattern='*.nc',full.names=TRUE)

##堆叠所有netCDF文件并提取变量“sst”

ncin_SST <- raster::stack(filenames, varname = "sst")

打印结果

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, ... 

这是一个 AQUA MODIS 文件的打印输出:

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: data@oceancolor.gsfc.nasa.gov
        publisher_email: data@oceancolor.gsfc.nasa.gov
        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

将 CRS 更改为 WG84/UTM34

ncin_SST_34 <- projectRaster(ncin_SST, crs = "+proj=utm + zone=34 + datum=WGS84")

错误信息

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

要打开包含栅格数据的 ncdf 文件,您只需执行以下操作

library(terra)
x <- rast("filename.nc")

如果每个文件有多个变量你可以select像这样

x <- rast("filename.nc", "SST")

(或者用 sds("filename.nc")

将它们全部打开为 SpatRasterDataset

如果你有一个向量 ff,其中包含同一区域(相同范围和分辨率)的多个 ncdf 文件的文件名,你可以这样做

x <- rast(ff, "SST")

虽然可以将您的数据转换为另一个坐标参考系统,例如像这样的 UTM

p <- project(x, "+proj=utm + zone=34 + datum=WGS84")

在大多数情况下这不是一个好主意。首先,您确实应该提供要投影到的模板栅格(请参阅 ?project),但更重要的是,投影栅格数据会扭曲值。所以在数据分析的时候尽量避免这种情况(大部分情况下做个图就可以了)。而是将空间矢量数据投影到栅格数据的 CRS。