时间戳

首页 > 解决方案 > 在 R 中导入 NetCDF 文件时，光栅砖中的经纬度丢失

r - 在 R 中导入 NetCDF 文件时，光栅砖中的经纬度丢失

问题描述

test <-brick(file.name) 我的数据非常大，很难上传。

我使用这段代码读入了一个 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

如果我使用 lat longs 拉出，我可以看到它们存在：

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

绘图以回应下面的评论

标签： rgisrasternetcdf

解决方案

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

该文件的坐标参考系 (crs) 不是 lon/lat。似乎它没有具体说明它是什么，或者，如果有，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)

推荐阅读