0

我有以下NetCDF文件 - 我正在尝试转换为栅格，但有些地方不对劲。没有给出 NetCDF 文件的投影，但基于我从它那里收到的软件，它应该是 LatLong，但可能是圆柱等面积。我尝试了两种方法，但我不断收到这种失真，这使得无法在正确的位置查询值。我知道网格的间距不均匀，不确定这是否会影响最终结果（这里是 ArcGIS 的视觉效果，但在 R 中这是同样的问题，除非使用 levelplot 函数绘制）。 在此处输入图像描述

library(raster)
library(ncdf4)
library(lattice)
library(RColorBrewer)

setwd("D:/Results")
climexncdf <- nc_open("ResultsSO_month.nc")

lon <- ncvar_get(climexncdf,"Longitude")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(climexncdf,"Latitude")
nlat <- dim(lat)
head(lat)

dname <- "Weekly Growth Index"

t <- ncvar_get(climexncdf,"Step")
tmp_array <- ncvar_get(climexncdf,dname)
tmp_stack <- vector("list",length(t))

for (i in 1:length(t)) {
        tmp_stack[[i]] <- tmp_array[,,i]
}

YearData <- vector("list",52)
for (i in 1:4) {
        YearData[[i]] <- tmp_array[,,i]
}   

Month1 <- YearData[c(1,2,3,4)]

# Calculate monthly averages
M1Avg <- Reduce("+",Month1)/length(Month1)

# Replace 0's with NA's
M1Avg[M1Avg==0] <- NA

# Piece of code that gives me what I need:
grid <- expand.grid(lon=lon, lat=lat)
cutpts <- seq(0,1,0.1)

# Convert to raster - work to include lat and long

M1Avg_reorder <- M1Avg[ ,order(lat) ]
M1Avg_reorder <- apply(t(M1Avg_reorder),2,rev)

M1AvgRaster <- raster(M1Avg_reorder,
                      xmn=min(lon),xmx=max(lon),
                      ymn=min(lat),ymx=max(lat),
                      crs=CRS("+proj=longlat +datum=WGS84"))
                        #crs=CRS("+proj=cea  +lat_0=0 +lon_0=0"))

r <- projectRaster(M1AvgRaster,crs=CRS("+proj=longlat +datum=WSG84"))

plot(M1AvgRaster)

# Location file not included but any locations can be entered
locations <- read.csv("Locations.csv", header=T)
coordinates(locations) <- c("y","x")

data <- extract(M1AvgRaster,locations)

writeRaster(M1AvgRaster, "M1AvgRaster_Globe_projWGSTest", format = "GTiff")
4

1 回答 1

1

python 版本显示，在重新排序后，至少数据的位置似乎是正确的。但是，数据文件看起来很奇怪，我看到 python netcdf 库中的数据实际上已经损坏，我以前从未见过有很多不同的 NetCDF 文件。此外，分块和压缩设置很奇怪，最好不要应用它们。

但是获取情节的最小python示例在这里：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset

ff = Dataset('ResultsSO_month.nc')
test_var = np.copy(ff.variables['Maximum Temperature'][:])
## reorder latitudes
latindex = np.argsort(ff.variables['Latitude'][:])
## Set up map and compute map coordinates
m = Basemap(projection='cea', llcrnrlat=-90, urcrnrlat=90, 
llcrnrlon=-180, urcrnrlon=180, resolution='c')
grid_coords = np.meshgrid(ff.variables['Longitude'[:],ff.variables['Latitude'][latindex])
X,Y = m(grid_coords[0],grid_coords[1])
## Plot
m.pcolormesh(X,Y,test_var[0,latindex,:])
m.drawcoastlines()
plt.colorbar()
plt.show()

在此处输入图像描述

于 2017-05-10T18:17:53.827 回答