我已经下载从NC文件
f=open.ncdf("file.nc")
[1] "file Lfile.nc has 2 dimensions:"
[1] "Longitude Size: 1440"
[1] "Latitude Size: 720"
[1] "------------------------"
[1] "file filr.nc has 8 variables:"
[1] "short ts[Latitude,Longitude] Longname:Skin Temperature (2mm) Missval:NA"
然后我想与变量soil_moisture_c工作
A = get.var.ncdf(nc=f,varid="soil_moisture_c",verbose=TRUE)
然后我绘制与image(A)
我得到了如下图所示的地图,我甚至调换它image(t(a))
但改为另一个方向,而不是它应该如何。 总之,要想知道什么是错的,我用的NetCDF观众甲胄和地图被正确地绘制成你可以看到下面。
正如补充@mdsumner远更好的解决办法,这里是一个办法做到这一点利用图书馆ncdf
只。
library(ncdf)
f <- open.ncdf("LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040101.nc")
A <- get.var.ncdf(nf,"soil_moisture_c")
所有你需要的是找到你的尺寸,以便有一个连贯的X和Y轴。 如果你看看你的netCDF对象的尺寸,在这里你所看到的:
str(f$dim)
List of 2
$ Longitude:List of 8
..$ name : chr "Longitude"
..$ len : int 1440
..$ unlim : logi FALSE
..$ id : int 1
..$ dimvarid : num 2
..$ units : chr "degrees_east"
..$ vals : num [1:1440(1d)] -180 -180 -179 -179 -179 ...
..$ create_dimvar: logi TRUE
..- attr(*, "class")= chr "dim.ncdf"
$ Latitude :List of 8
..$ name : chr "Latitude"
..$ len : int 720
..$ unlim : logi FALSE
..$ id : int 2
..$ dimvarid : num 1
..$ units : chr "degrees_north"
..$ vals : num [1:720(1d)] 89.9 89.6 89.4 89.1 88.9 ...
..$ create_dimvar: logi TRUE
..- attr(*, "class")= chr "dim.ncdf"
因此,你的尺寸为:
f$dim$Longitude$vals -> Longitude
f$dim$Latitude$vals -> Latitude
现在,您Latitude
90反潜至-90相反,后者的image
更喜欢,所以:
Latitude <- rev(Latitude)
A <- A[nrow(A):1,]
最后,你发现,你的对象A的X和Y翻转,所以你需要转它,你的NA值表示由值某些原因-32767
:
A[A==-32767] <- NA
A <- t(A)
最后的情节:
image(Longitude, Latitude, A)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)
编辑:要做到这一点在你31页的文件,让我们打电话给你的文件名的矢量ncfiles
和yourpath
您保存它们的目录(为简单起见我要去假设你的变量总是叫soil_moisture_c
和您的NAS总是-32767
):
ncfiles
[1] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040101.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040102.nc"
[3] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040103.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040104.nc"
[5] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040105.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040106.nc"
[7] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040107.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040108.nc"
[9] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040109.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040110.nc"
[11] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040111.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040112.nc"
[13] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040113.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040114.nc"
[15] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040115.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040116.nc"
[17] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040117.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040118.nc"
[19] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040119.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040120.nc"
[21] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040121.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040122.nc"
[23] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040123.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040124.nc"
[25] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040125.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040126.nc"
[27] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040127.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040128.nc"
[29] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040129.nc" "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040130.nc"
[31] "LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T173800Z_20040131.nc"
yourpath
[1] "C:\\Users"
library(ncdf)
library(maptools)
data(wrld_simpl)
for(i in 1:length(ncfiles)){
f <- open.ncdf(paste(yourpath,ncfiles[i], sep="\\"))
A <- get.var.ncdf(f,"soil_moisture_c")
f$dim$Longitude$vals -> Longitude
f$dim$Latitude$vals -> Latitude
Latitude <- rev(Latitude)
A <- A[nrow(A):1,]
A[A==-32767] <- NA
A <- t(A)
close.ncdf(f) # this is the important part
png(paste(gsub("\\.nc","",ncfiles[i]), "\\.png", sep="")) # or any other device such as pdf, jpg...
image(Longitude, Latitude, A)
plot(wrld_simpl, add = TRUE)
dev.off()
}
原因是,您所使用的的NetCDF接口是非常低的水平,并已完成所有被读出的变量没有它的任何维度信息。 网格的方向实在是任意的,并且需要在特定的上下文中理解的坐标信息。
library(raster) ## requires ncdf package for this file
d <- raster("LPRM-AMSR_E_L3_D_SOILM3_V002-20120520T185959Z_20040114.nc", varname = "soil_moisture_c")
(我用不同的文件到你的,但它应该工作一样)。
事实证明,即使栅格没有得到工作没有这个权利,但它可以很容易地纠正:
d <- flip(t(d), direction = "x")
这换位数据和翻转围绕“X”(经度),保持地理坐标从原来的环境。
从maptools检查地图绘制起来:
plot(d)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)
还有许多其他的方式通过读取文件中的维度信息来实现这一点,但是这至少是一个快捷方式做最艰难的工作适合你。
你也可以简单地颠倒使用CDO命令行纬度:
cdo invertlat file.nc file_inverted.nc