如何将我的时空 NetCDF 数据转换为空间数据?

How to convert my spatiotemporal NetCDF data to spatial data?

我是 R 的初学者,我完全被这个问题困扰了。你可以从下面的link下载netCDF文件看看。

https://drive.google.com/file/d/0ByY3OAw62EShbkF6VWNFUkRYMmM/view?usp=sharing

^这是我的NetCDF大气数据文件,有8个变量和8个维度。在这里,我感兴趣的变量是:

TIMSID为站点编号(包括城市站点、农村站点等)
URBAN :: 城市站点数 [urban 是 3 行 250 列的矩阵。第1行是城市站点数,第2行是纬度,第3行是经度。]
TIME :: 数据收集于 2012 年 3 月 1 日至 2012 年 5 月 ['time' 的编码为 YYYYMMDDHH]
PM10 :: 每个地点每个站测量的每小时颗粒物浓度

我只需要处理这个大数据集中的这 4 个变量。

我只需要分离“2012 年 3 月 1 日”城市站点的 PM10 值数据。 (实际上我需要在 TIMSID 变量中找到哪些站点是城市站点,并匹配 2012 年 3 月 1 日城市站点的相应 PM10 值。)

例如,在 TIMSID 中,存在城市、农村等不同类型的站点,名称为 111121、111122、111123、111124,但城市站点编号为 111121、111123..等,因此我必须仅考虑来自 TIMSID 的城市站点data 并且要匹配对应的pm10值,时间,纬度,经度。最后想制作一个新的数据集。

最后的table/dataset应该是~column1-time(only 1 march 2012),column2-urban sites number, Column (3,4)-对应城市站点的经纬度,Column 5- hourly各城市站点的pm10值

我已经使用以下这些命令从 NetCDF 文件中读取数据。但是我不明白我应该做些什么...

install.packages("ncdf",dependencies=TRUE)
library(ncdf)

nc<-open.ncdf("2012_03_05_PM10_surface.nc")
print(nc)

tmsid<-get.var.ncdf(nc,"TMSID")
timsid

urban<-get.var.ncdf(nc,"urban")
urban
time<-get.var.ncdf(nc,"TIME")

pm10<-get.var.ncdf(nc,"PM10")

由于我是R的初学者,所以我只知道基本的命令。我想不通,我应该学习哪个特定的包来解决这个问题。请帮帮我好吗?提前感谢您的宝贵时间。如果您需要任何进一步的信息,请随时问我。

library(ncdf)
nc <- open.ncdf("2012_03_05_PM10_surface.nc")
tmsid <- get.var.ncdf(nc,"TMSID")
urban <- get.var.ncdf(nc,"urban")
time <- get.var.ncdf(nc,"TIME")
pm10 <- get.var.ncdf(nc,"PM10")

先来看看nc:

[1] "file ~/Downloads/2012_03_05_PM10_surface.nc has 8 dimensions:"
[1] "data_num   Size: 683016"
[1] "ncl1   Size: 683016"
[1] "obsnum_urban   Size: 250"
[1] "ID_LAT_LON   Size: 3"
[1] "obsnum_road   Size: 33"
[1] "obsnum_background   Size: 5"
[1] "obsnum_rural   Size: 16"
[1] "ncl7   Size: 683016"
[1] "------------------------"
[1] "file ~/Downloads/2012_03_05_PM10_surface.nc has 8 variables:"
[1] "int TMSID[data_num]  Longname:TMSID Missval:NA"
[1] "int TIME[ncl1]  Longname:TIME Missval:NA"
[1] "float PM10[data_num]  Longname:PM10 Missval:1e+30"
[1] "float urban[ID_LAT_LON,obsnum_urban]  Longname:urban Missval:1e+30"
[1] "float road[ID_LAT_LON,obsnum_road]  Longname:road Missval:1e+30"
[1] "float background[ID_LAT_LON,obsnum_background]  Longname:background Missval:1e+30"
[1] "float rural[ID_LAT_LON,obsnum_rural]  Longname:rural Missval:1e+30"
[1] "int TMS_JULIAN[ncl7]  Longname:TMS_JULIAN Missval:NA"

它告诉我们的是urban的行是ID,纬度和经度。然后我们有 tmsid 给出与 time 向量大小相同的 ID 向量:每个 data_num,i。 e. PM10 中每个数据点的一对 ID 时间,这意味着我们将能够通过 ID(由 urban 的第一行给出)和时间戳(从 2012030101 到2012030124).

# First we need to make a dataframe out of urban, for convenience.
urban <- as.data.frame(t(urban))
colnames(urban) <- c("ID", "LAT", "LON")
# Then we do the subsetting using a lapply, so we can batch-subset:
res <- lapply(urban$ID, 
              function(x)data.frame(ID=x,
                                    pm=pm10[tmsid%in%x & time%in%2012030101:2012030124], 
                                    time=2012030101:2012030124))
# Which gives us a list of sub-dataframes that we want to compress back into a single dataframe:
res <- do.call(rbind,res)
# Finally we merge that with the original urban dataframe
# so that each entry has its own LAT and LON:
res <- merge(res, urban, by="ID")
res
#         ID   pm       time      LAT      LON
#1    111121   42 2012030101 37.56464 126.9760
#2    111121   36 2012030102 37.56464 126.9760
#3    111121   46 2012030103 37.56464 126.9760
#4    111121   40 2012030104 37.56464 126.9760
#5    111121   36 2012030105 37.56464 126.9760
#...
#5995 831154   81 2012030119 37.52662 126.8064
#5996 831154   72 2012030120 37.52662 126.8064
#5997 831154   81 2012030121 37.52662 126.8064
#5998 831154   70 2012030122 37.52662 126.8064
#5999 831154   74 2012030123 37.52662 126.8064
#6000 831154   74 2012030124 37.52662 126.8064

250 个城市站点 X 24 小时 = 6 000 个数据点,这确实是我们在这里得到的。