GPM 降雨量数据处理 -R(坐标系转换)
背景
今天给大家介绍下,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)
image.png
参考
1.Geocomputation with R 2.Experiments with sf (spatial data in r) 3.Rasterizing an sf vector object
相关文章
- JS转换HTML转义符,防止javascript注入攻击,亲测可用「建议收藏」
- C#*.bmp,*.jpg,*.png指定颜色,转换为透明背景的png[通俗易懂]
- dotnet 将任意时区的 DateTimeOffset 转换为中国时区时间文本
- Java 把 Map 的值(Value)转换为 Array, List 或 Set
- Python的时间转换详解编程语言
- [javaSE] 进制转换(二进制十进制十六进制八进制)详解编程语言
- 甲骨文云计算平台推出免费迁移转换工具 吸引企业从微软或亚马逊迁移过来
- Oracle数据转换成字符—实现快速而可靠的数据处理(oracle转换为字符)
- 轻松转换格式,Linux转码工具下载推荐(linux转码工具下载)
- MySQL中如何进行IP地址转换(mysql中ip转换)
- dtcms成功转换为oracle(dtcms oracle)
- csv数据快速转换至MySQL数据库(.csv转mysql)
- c#转换全角半角方法示例
- python转换摩斯密码示例