美文网首页遥感
GPM 降雨量数据处理 -R(坐标系转换)

GPM 降雨量数据处理 -R(坐标系转换)

作者: jamesjin63 | 来源:发表于2020-06-11 22:00 被阅读0次

背景

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

这里需要强调的一点就是,降雨数据主要在NASA网站主要包括TRMMGPM项目
下载的数据是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

相关文章

  • GPM 降雨量数据处理 -R(坐标系转换)

    背景 今天给大家介绍下,R处理NASA下载的降雨量数据在进行环境数据分析时候,经常需要用到降雨量的信息,而NASA...

  • 2018-03-04

    常用坐标系统知识点 1.坐标系统之间的转换 (1)坐标系分类 不同参心坐标系之间的转换、不同地心坐标系之间的转换;...

  • C#:霍夫变换圆检测

    一、霍夫圆检测原理: 与霍夫直线检测相似,从平面坐标系到极坐标系转换三个参数 C( x0,y0,r ) 只是参数方...

  • 地图坐标转换

    地图坐标转换 简介 各地图API坐标系统比较与转换; WGS84坐标系:即地球坐标系,国际上通用的坐标系。设备一般...

  • 高等数学预备知识

    极坐标 极坐标系是由极轴、极径组成极坐标系上的点表示为(ρ,θ)极坐标系上的点转换成直角坐标直角坐标系上的点转换成...

  • week51 坐标变换与坐标系变换

    坐标转换是一个坐标在不同坐标系下的表示,而坐标系转换不同坐标系的相对位姿关系。 TF介绍TF(TransForm)...

  • 地理坐标系转换 API 接口

    地理坐标系转换 API 接口 提供地理信息坐标系的相互转换。 1. 产品功能 支持多种地理信息坐标系; 高精度坐标...

  • 顶视图-世界坐标系-相机坐标系-Scara坐标系转换

    顶视图-世界坐标系-相机坐标系-Scara坐标系转换 1. 定义 Scara相机坐标系 标准相机坐标系 顶视图坐标...

  • Qt 绘图转换

    转换 QTransform 用于指定坐标系的 2D 转换 - 平移、缩放、扭曲(剪切)、旋转或投影坐标系。绘制图形...

  • 环境

    1. 导入需要的库 2. 定义三角函数用于坐标系转换 3. 定义转换矩阵:大地坐标系---->随体坐标系 4. 定...

网友评论

    本文标题:GPM 降雨量数据处理 -R(坐标系转换)

    本文链接:https://www.haomeiwen.com/subject/bzbitktx.html