如何将时空NetCDF数据转换成空间数据？

Question

我是R.的初学者，我对这个问题完全着迷。您可以从下面的链接下载netCDF文件来查看。

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

这是我的NetCDF大气数据文件，包含8个变量和8个维度。在这里，我感兴趣的变量是：

TIMSID是指站点的数量(包括城市站点、农村站点等)。

城市:城市用地数量为3列250列矩阵。row1是城市站点的数量，row2是纬度，第3行是经度。

时间:数据收集于2012年3月1日至2012年5月，编码时间为YYYYMMDDHH

PM10 :：在每个站点的每个站测量每小时颗粒物浓度

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

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

例如，在TIMSID中，城市、农村等不同类型的站点被命名为111121,111122,111123,111124，但城市站点号是111121,111123..etc，所以我必须只考虑来自TIMSID数据的城市站点，并且需要匹配相应的pm10值、时间、纬度、经度。最后想要做一个新的数据集。

最后的表/数据集应该是~列1-时间(2012年3月1日)，第2列-城市场址数、列(3,4)-latitude和相应城市站点的经度，每个城市站点的每5小时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的初学者，所以我只知道基本的命令。我不知道，我应该学会解决这个问题的具体方案。请帮帮我?提前感谢您宝贵的时间。如果您需要进一步的信息，请随时问我。

Answer 1

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向量:1个data_num，即PM10中的每个数据点一对ID时间，这意味着我们可以根据ID(由urban的第一行给出)和时间戳(从2012030101到2012030124)对pm10进行子集。

# 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个数据点，这就是我们在这里得到的。