-1

NetCDFファイルで逆座標変換を行うにはどうすればよいですか？75の経度値と36の緯度値を持つグリッドがあります。

nc<-create.n("filename.nc")
#Dimentions
dim.def.nc(nc,"lon",75)
dim.def.nc(nc,"lat",36)
dim.def.nc(nc,"time",365)
#Vars
var.def.nc(nc,"Observation","NC_FLOAT", c(1,0,2))
var.def.nc(nc,"lon","NC_FLOAT", c(0))
var.def.nc(nc,"lat","NC_FLOAT", c(1))
var.def.nc(nc,"time","NC_FLOAT", c(2))
(...)

unidataのドキュメントによると、netcdfで（lat、lon）から（x、y）への逆変換を実行できるはずですが、これをどのように実行できるかわかりません。緯度経度グリッドをランベルト正角円錐図法グリッドに変換したいと思います。

4

1 に答える 1

4

これは、netCDFファイルを再投影する方法です。基本的に、必要な経度、緯度、およびデータを取得し、パッケージrgdalで使用できるシェープファイルを作成spしますmaptools。（ここでの例では、ここでダウンロードしたNOAAのデータを使用しています）

library(ncdf)
library(rgdal)
library(sp)
library(maptools)
nc <- open.ncdf("20130128-ABOM-L4HRfnd-AUS-v01-fv01_0-RAMSSA_09km.nc")

# Grab the longitude, latitude and data
lon <- nc$dim$lon$vals
lat <- nc$dim$lat$vals
sst <- get.var.ncdf(nc,"analysed_sst")

# Create a SpatialPointsDataFrame object
lonlat <- expand.grid(lon,lat)
sst <- as.data.frame(matrix(sst,ncol=1))
dat <- SpatialPointsDataFrame(lonlat, data=sst,
                              proj4string=CRS("+proj=longlat +datum=WGS84 "))

# And then reproject
dat2 <- spTransform(dat,CRS("+proj=lcc"))  
# Of course you have to write the proj4 string that corresponds exactly to the desired projection.
于 2013-01-29T16:44:36.160 に答える