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 に答える