用R处理NASA数据（.hdf 或.nc文件）

Jamesjin63

发布于 2022-10-25 14:41:31

1.2K0

发布于 2022-10-25 14:41:31

文章被收录于专栏：EpiHub

1.下载NASA数据

这里不在赘述，参考如何获取NASA数据，下面的例子根据下载的LandCover与Rainfall数据进行展示，如何利用R语音进行读取，然后绘图。先加载所需R包及地图文件

library(ncdf4)
library(rgdal)
library(gdalUtils)
library(raster)
library(rasterVis)
library(sf)
library(exactextractr)
library(maptools)
library(cleangeo)
library(rworldmap)
rm(list=ls())
# continental shapefile
cont = getMap()
cont = clgeo_Clean(cont)
cont = sapply(levels(cont$continent),
               FUN = function(i){
                 poly = gUnionCascaded(subset(cont, continent == i))
                 poly = spChFIDs(poly, i)
                 SpatialPolygonsDataFrame(poly, data.frame(continent = i, row.names = i))
               }, USE.NAMES = TRUE)
cont = Reduce(spRbind, cont)
# read shp files
setwd("/Users/Desktop/NASA/LandCover")
CHN = st_read("/Users/Desktop/NASA/LandCover/shp/最新的全国区划/县.shp")
CHN_sp = readOGR("/Users/Desktop/NASA/LandCovershp/最新的全国区划/县.shp")

2.读取hdf文件

将hdf文件存在Landcover文件夹目录下，然后查看hdf文件

> hdf=list.files(pattern = ".hdf")
> hdf
[1] "MCD12C1.A2010001.006.2018053185051.hdf" "MCD12C1.A2011001.006.2018053185321.hdf"
[3] "MCD12C1.A2014001.006.2018053185556.hdf" "MCD12C1.A2015001.006.2018053185652.hdf"
[5] "MCD12C1.A2016001.006.2018324172410.hdf" "MCD12C1.A2017001.006.2019192025407.hdf"
[7] "MCD12C1.A2018001.006.2019200161458.hdf"

譬如我们需要读取第二个文件 MCD12C1.A2011001.006.2018053185321.hdf；这里的gdalinfo(hdf_filesname) 可以显示该hdf文件的详细列表信息，经纬度，坐标系，年份及卫星信息；sds就是我们想要的数据，其中Majority_Land_Cover_Type_1是根据MCD12C1第一个分类标准，将地球的植被覆盖分成25个类型；具体见官网说明文档。

hdf_filesname=hdf[2]
hdf_tif_name=paste0(unlist(str_split(hdf_filesname,".hdf"))[1],".tif")
gdalinfo(hdf_filesname)
hdf_time = str_extract(hdf_filesname,"()[0-9]{7}") 
sds = get_subdatasets(hdf_filesname)
> sds
[1] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_1"           
[2] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_1_Assessment"
[3] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Land_Cover_Type_1_Percent"            
[4] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_2"           
[5] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_2_Assessment"
[6] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Land_Cover_Type_2_Percent"            
[7] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_3"           
[8] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Majority_Land_Cover_Type_3_Assessment"
[9] "HDF4_EOS:EOS_GRID:MCD12C1.A2011001.006.2018053185321.hdf:MOD12C1:Land_Cover_Type_3_Percent"

最关键的一步，就是读取的第一个Majority_Land_Cover_Type_1文件，从hdf抽取出来转换成tiff文件。你会发现，你的文件夹下多了个相同hdf名字的tiff文件。然后读取tiff到raster就可以了

gdal_translate(sds[1], dst_dataset = hdf_tif_name)  # change hdf to tiff
hdf_raster=raster(hdf_tif_name)                     # read tiff as raster
# covert F into T
names(hdf_raster)=hdf_time
hdf_raster
class      : RasterLayer 
dimensions : 3600, 7200, 25920000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /Users/Desktop/NASA/Landcover/MCD12C1.A2011001.006.2018053185321.tif 
names      : X2011001 
values     : 0, 255  (min, max)

3.绘图

hdf_raster是我们提取到的25类 landcover，接下来就是绘图部分

rasterVis::levelplot(hdf_raster, margin = NA, par.settings = RdBuTheme) + 
  layer(sp.polygons(cont))

屏幕快照 2020-06-03 下午3.55.26.png

4.提取中国范围内的数据

接下来就是根据中国地图，将Landcover裁剪至China map

# crop within China
CHN_cropped = crop(hdf_raster, CHN)
CHN_masked = mask(hdf_raster, CHN) # take >5mins
# plot
mapTheme = rasterTheme(region=rev(brewer.pal(8,"Greens")))
rasterVis::levelplot(CHN_cropped, margin = NA, par.settings = mapTheme) +
  layer(sp.polygons(CHN_sp1))

屏幕快照 2020-06-03 下午3.59.05.png

接下来，我们用ggplot展示下结果。(制图反应时间较长)

第一种方法，加载SpatialPolygonsDataFram地图

第二种方法，加载Classes ‘sf’格式地图

#ggplot with raster
# change raster into dataframe
df_CHN_masked=as.data.frame(CHN_masked,xy=T) 
colnames(df_CHN_masked)=c("x","y","LandCover")
  
# change SpatialPolygonsDataFram into dataframe
CHNshp = fortify(CHN_sp)

# Method 1
ggplot(df_CHN_masked  %>% na.omit() ) +
  geom_raster(aes(x,y,fill=LandCover))+
  scale_fill_gradientn(colours=c("grey","green")) +
  coord_quickmap()+
  geom_polygon(data = CHNshp, aes(long, lat, group = group),
             fill=NA,color="grey50", size=0.1)

# Method 2
ggplot(df_CHN_masked  %>% na.omit() ) +
  geom_tile(aes(x,y,fill=LandCover))+
  scale_fill_gradientn(colours=c("grey","blue")) +
  geom_sf(data=CHN,size=1,fill=NA) 

# with the county seat
#source1
Chinamap=read_sf("https://gw.alipayobjects.com/os/rmsportal/JToMOWvicvJOISZFCkEI.json")
ggplot()+
  geom_sf(data=Chinamap)

ggplot(df_CHN_masked  %>% na.omit() ) +
  geom_tile(aes(x,y,fill=LandCover))+
  scale_fill_gradientn(colours=c("grey","blue")) +
  geom_sf(data=Chinamap,size=0.2,fill=NA,color="black")