######################################################            

#  This script is used R to read netcdf data        

#           and export to asc            

#              Author: Mcr                                  

#nc多层数据没能在GIS中打开，将nc转存为多个asc后可以使用。

rm(list=ls())

library(lattice)

setwd("...") #设定nc数据所在位置路径

print.nc(my.data)               #查看内部属性

#n.depth

r <- raster(ncol=7560, nrow=4320,xmn=73, xmx=136, ymn=18, ymx=54)

#建立一个空栅格图，7560列，4320行，经度73-136，纬度18-54

TH33_01_tmp<-var.get.nc(my.data,"TH33",c(1,1,1),c(7560,4320,1))#读取第一层文件

#☆c(1,1,1) 表示c（读取nc栅格图起始单元格的行位置，读取nc栅格图起始单元格的列位置，读取nc栅格图

#起始单元格的层位置）

#☆c(7560,4320,1) 表示c(读取nc栅格图单元格总行数，读取nc栅格图单元格总列数，读取nc栅格图层数）

TH33_01_tmp<-TH33_01_tmp[,4320:1] #将列顺序倒置，这与输出图的方向有关。

plot(r,asp=1) #展示r， 注意图的方向

TH33_01<-as.data.frame(r) #将r栅格内数据转为data.frame

TH33_01[is.na(TH33_01)]=-9999 #将NA替换为-9999，要不然asc数据读入到GIS会出错。

writeLines('ncolstt7560nnrowstt4320tnxllcornert73.00417nyllcornert18.00417ncellsizet0.00833nNODATA_valuet-9999', fileCon)

#为要输出的asc文件写上表头，包括列数，行数，起始经度纬度，分辨率和空值数据用-9999表示。

write.table(TH33_01,'TH33_01.asc',append=TRUE,sep=" ",row.names=FALSE,col.names=FALSE)

#文中代码只输出了七层数据得第一层数据，第i层用c(1,1,i)，c(7560,4320,1)。

#输出范围可以通过行列数进行调整，而相应的经纬度也需要调整。

###文中会有纰漏，请见谅指正。