背景

今天给大家介绍下，R处理NASA下载的降雨量数据
在进行环境数据分析时候，经常需要用到降雨量的信息，而NASA提供了每年，每个月甚至每天的降雨量数据。
如何下载NASA降雨量数据，见此链接。

这里需要强调的一点就是，降雨数据主要在NASA网站主要包括TRMM与GPM项目
下载的数据是HDF5格式，如何在R读取HDF与tc文件，戳here。
TRAMM与GRM下载的HDF5格式在R中，会出现坐标与我们常用坐标系不一致的情况，
主要投影坐标系不同。

所以这篇文章，这要介绍raster如何转换成常规的4236坐标系。

1.读取数据

library(ncdf4)
library(rgdal)
library(gdalUtils)
library(raster)
library(rasterVis)
library(sf)
library(exactextractr)
library(maptools)
library(cleangeo)
library(rworldmap)
rm(list=ls())
# read shp files
# 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)

# set to GPMM directory
file.remove(list.files(pattern = ".tif"))
hdf=list.files(pattern = ".HDF5")

# read HDF5
hdf_filesname=hdf[1] #first one
gdalinfo(hdf_filesname) 
sds = get_subdatasets(hdf_filesname)
sds

# select varibles: precipitation
gdal_translate(sds[4], dst_dataset = hdf_tif_name)# change hdf to tiff
#read tiff as raster
hdf_raster=raster(hdf_tif_name)

上述主要是将HDF5文件转换成Raster文件，找到储存在HDF5文件中的precipitation位置。然后存储到hdf_raster当中。

2.Raster转换

接下来是关键性的一步，过程比较长。cont是世界地图的SpatialPolygonsDataFrame 数据，我们在前面加载好
我们先看看hdf_raster长什么样子。

image.png

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

image.png

嚯嚯，这里的hdf_raster与左下角的cont一点也不对应，怎么办？
我们将hdf_raster旋转一下，这样子可以看到差不多正常了。
但是cont还是在左下角，坐标对应不上。

# rotate the x and y axis
s2 = t(flip(hdf_raster, direction='y' ))
# plot
rasterVis::levelplot(s2, margin = NA, par.settings = RdBuTheme) + layer(sp.polygons(cont))

image.png

下面我们新建一个对应lat与long的空的Raster为rasterNoProj，可以指定分辨率为0.5.
将rasterNoProj 转换成数据库data.frame，包含了x，y坐标信息。
然后我们之前旋转后的s2也转换data.frame格式。

#craet new raster with 360*720,0.5res
newMatrix  = matrix(runif(3600*1800,100,999), nrow = 1800, ncol =3600)
rasterNoProj = raster(newMatrix,xmn = -180, xmx = 180, ymn = -90, ymx = 90)

# change the rasterNoProj x,y to TRMM_raster x,y
spg = as.data.frame(rasterNoProj, xy=TRUE)
TRMM_spg = as.data.frame(s2, xy=TRUE)

接下来就是TRMM_spg 数据放在spg数据框里面。

# Correct extent
spg$layer=TRMM_spg$layer
coordinates(spg) = ~ x+y
gridded(spg) = TRUE
rasterDF = raster(spg);rm(spg,newMatrix,TRMM_spg);rm(rasterNoProj)
# crs(rasterDF)="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 "

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

image.png

3. sf方法(耗时太长，不建议尝试)

其实sf方法更简单。但是s2数据太大，转换成sf时间较长，
先喝口水。慢慢等待。
缺点，在制图过程中，也需要很长时间才能出图。

# rgdal::make_EPSG() %>% 
#   DT::datatable()
# change to sf
df_sf =s2 %>% rasterToPoints(spatial = T) %>%
  st_as_sf()
# change to 4326
st_crs(df_sf) = 4326

# plot
cont_sf=st_as_sf(cont)
ggplot() +
  geom_sf(data=cont_sf,size=0.2,fill=NA)+
  geom_sf(data=df_sf)