在R中导入NetCDF文件时丢失了光栅砖的长整型

Question

测试

我使用此代码读入NetCDF文件。 NetCDF文件包含文件中的时间较长。我知道这一点是因为我可以提取它们（见下文）。但是，当我以砖块形式读取netCDF文件时，砖块的范围不能反映netCDF文件中的经时长度。它似乎只是将栅格居中于美国大陆的中心（请参见下面的第一幅图）：

library(ncdf4)
library(raster)

test <- brick(file.name)
test

class      : RasterBrick 
dimensions : 141, 240, 33840, 7273  (nrow, ncol, ncell, nlayers)
resolution : 25000, 25000  (x, y)
extent     : -3e+06, 3e+06, -1762500, 1762500  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : C:/Users/anfangen/Box/Environmental-Data/Model-Pollen-Data/Raw-Data/RAGcount_1997-2016.nc 
names      : X1997.01.02.00.00.00, X1997.01.03.00.00.00, X1997.01.04.00.00.00, X1997.01.05.00.00.00, X1997.01.06.00.00.00, X1997.01.07.00.00.00, X1997.01.08.00.00.00, X1997.01.09.00.00.00, X1997.01.10.00.00.00, X1997.01.11.00.00.00, X1997.01.12.00.00.00, X1997.01.13.00.00.00, X1997.01.14.00.00.00, X1997.01.15.00.00.00, X1997.01.16.00.00.00, ... 
Date/time  : 1997-01-02 00:00:00, 2016-12-29 00:00:00 (min, max)
varname    : count 

如果我使用来拉出经纬度，我可以看到它们存在：

test_nc_open <- nc_open(file_name.nc')

## Get the lat longs for the ragweed file. These are in matrix form. It is as yet unclear to me why we need them in this form. 
lat <- ncvar_get(test_nc_open, varid = 'xlat')
lon <- ncvar_get(test_nc_open, varid = 'xlon')

lat

           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]    [,10]    [,11]    [,12]
  [1,] 18.86005 19.05895 19.25808 19.45742 19.65697 19.85674 20.05672 20.25690 20.45730 20.65790 20.85870 21.05970
  [2,] 18.92987 19.12906 19.32847 19.52810 19.72794 19.92799 20.12826 20.32874 20.52942 20.73031 20.93141 21.13270
  [3,] 18.99921 19.19868 19.39837 19.59828 19.79841 19.99875 20.19931 20.40007 20.60105 20.80223 21.00361 21.20520
  [4,] 19.06806 19.26781 19.46778 19.66797 19.86839 20.06901 20.26985 20.47091 20.67217 20.87363 21.07531 21.27719
  [5,] 19.13641 19.33644 19.53669 19.73717 19.93786 20.13877 20.33989 20.54123 20.74278 20.94453 21.14649 21.34866
  [6,] 19.20427 19.40458 19.60511 19.80586 20.00683 20.20802 20.40943 20.61105 20.81288 21.01492 21.21716 21.41962
  [7,] 19.27163 19.47221 19.67302 19.87405 20.07530 20.27677 20.47846 20.68036 20.88247 21.08479 21.28732 21.49006

我试图手动设置范围，但是当我绘制结果时，它与其他空间层不匹配：

## Get state boundaries 
states <- shapefile("GIS Data/State_boundaries/cb_2018_us_state_20m.shp")

test <- brick(file.name)
[1] "vobjtovarid4: error #F: I could not find the requsted var (or dimvar) in the file!"
[1] "var (or dimvar) name: crs"
[1] "file name: file.name.nc"


test <-
  setExtent(test, c(min(apply(lon, 2, min)), 
                    max(apply(lon, 2, max)), 
                    min(apply(lat, 2, min)), 
                    max(apply(lat, 2, max)) ), keepres =
              FALSE, snap = FALSE)

crs(test) <- crs(states)

plot(test[[1]])
plot(states, add=TRUE)

我不明白我在做什么错。我已经看过其他几个类似的问题，无法使用这些信息来解决这个问题。

print(test)

     3 variables (excluding dimension variables):
        float xlat[jx,iy]   
            long_name: Latitude on Cross Points
            standard_name: latitude
            units: degrees_north
            grid_mapping: crs
        float xlon[jx,iy]   
            long_name: Longitude on Cross Points
            standard_name: longitude
            units: degrees_east
            grid_mapping: crs
        float count[jx,iy,time]   
            cell_methods: time: point
            grid_mapping: crs
            coordinates: xlat xlon
            units: particles m-3
            standard_name: near-surface_particle_count
            long_name: Near-surface particle count

     3 dimensions:
        time  Size:7273   *** is unlimited ***
            long_name: time
            standard_name: time
            units: hours since 1949-12-01 00:00:00 UTC
            calendar: gregorian
            bounds: time_bnds
        iy  Size:141
            long_name: y-coordinate in Cartesian system
            standard_name: projection_y_coordinate
            units: m
            axis: Y
            _CoordinateAxisType: GeoY
        jx  Size:240
            long_name: x-coordinate in Cartesian system
            standard_name: projection_x_coordinate
            units: m
            axis: X
            _CoordinateAxisType: GeoX

Answer 1

尚不清楚您为什么要做自己的事情。在这种情况下，您要寻找的解决方案取决于您的目标。

该文件具有非lon / lat的坐标参考系（crs）。似乎它没有指定它是什么，或者，如果是，则raster无法理解它。无论您是谁提供文件的，都可能获得crs。要查看文件中是否有信息，可以执行print(test)

文件还存储单元格的经度和纬度。但是这些不在常规网格上，因此您不能使用它们。如果您确实需要栅格数据为lon / lat，则可以使用

x <- projectRaster(test, "+proj=longlat +datum=WGS84")

但是这会导致精度损失，因此最好避免。要从栅格数据中按点或面提取值，应首先将其转换为栅格数据的crs。您可以使用sp::spTransform。之后，您可以使用raster::extract

因此，使用SpatialPoints p，您可以执行类似的操作（但是使用正确的crs，在这种情况下utm 43肯定是错误的）

 putm <- spTransform(p, CRS("+proj=utm +zone=43 +datum=WGS84")) 
 e <- extract(test, putm)