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")
屏幕快照 2020-06-03 下午8.52.39.png
网友评论