在不规则网格上绘制数据的有效方法
stm*_*4tt 10 r netcdf ggplot2 map-projections r-sf
我使用组织在不规则二维网格上的卫星数据,其尺寸为扫描线(沿轨道尺寸)和地面像素(跨越轨道尺寸).每个中心像素的纬度和经度信息存储在辅助坐标变量中,以及四个角坐标对(纬度和经度坐标在WGS84参考椭球上给出).数据存储在netCDF4文件中.
我想要做的是在投影地图上有效地绘制这些文件(可能还有文件的组合 - 下一步!).
我的做法,到目前为止,灵感来自杰里米的Voisey的回答这个问题,一直在打造我的兴趣可变连杆的像素边界的数据帧,并使用ggplot2与geom_polygon用于实际的情节.
让我来说明我的工作流程,并提前为天真的方法道歉:我刚开始用一周或两周的R编码.
注意
要完全重现问题:
1.下载两个数据帧:so2df.Rda(22M)和pixel_corners.Rda(26M)
2.在您的环境中加载它们,例如
so2df <- readRDS(file="so2df.Rda")
pixel_corners <- readRDS(file="pixel_corners.Rda")
- 跳转到"合并数据帧"步骤.
初始设置
我要从我的文件中读取数据和纬度/经度边界.
library(ncdf4)
library(ggplot2)
library(ggmap)
# set path and filename
ncpath <- "/Users/stefano/src/s5p/products/e1dataset/L2__SO2/"
ncname <- "S5P_OFFL_L2__SO2____20171128T234133_20171129T003956_00661_01_022943_00000000T000000"
ncfname <- paste(ncpath, ncname, ".nc", sep="")
nc <- nc_open(ncfname)
# save fill value and multiplication factors
mfactor = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column",
"multiplication_factor_to_convert_to_DU")
fillvalue = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column",
"_FillValue")
# read the SO2 total column variable
so2tc <- ncvar_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column")
# read lat/lon of centre pixels
lat <- ncvar_get(nc, "PRODUCT/latitude")
lon <- ncvar_get(nc, "PRODUCT/longitude")
# read latitude and longitude bounds
lat_bounds <- ncvar_get(nc, "GEOLOCATIONS/latitude_bounds")
lon_bounds <- ncvar_get(nc, "GEOLOCATIONS/longitude_bounds")
# close the file
nc_close(nc)
dim(so2tc)
## [1] 450 3244
因此对于这个文件/传递,3244扫描线中的每一个都有450个地面像素.
创建数据帧
在这里,我创建了两个数据帧,一个用于值,一些用于后处理,一个用于纬度/经度边界,然后合并两个数据帧.
so2df <- data.frame(lat=as.vector(lat), lon=as.vector(lon), so2tc=as.vector(so2tc))
# add id for each pixel
so2df$id <- row.names(so2df)
# convert to DU
so2df$so2tc <- so2df$so2tc*as.numeric(mfactor$value)
# replace fill values with NA
so2df$so2tc[so2df$so2tc == fillvalue] <- NA
saveRDS(so2df, file="so2df.Rda")
summary(so2df)
## lat lon so2tc id
## Min. :-89.97 Min. :-180.00 Min. :-821.33 Length:1459800
## 1st Qu.:-62.29 1st Qu.:-163.30 1st Qu.: -0.48 Class :character
## Median :-19.86 Median :-150.46 Median : -0.08 Mode :character
## Mean :-13.87 Mean : -90.72 Mean : -1.43
## 3rd Qu.: 31.26 3rd Qu.: -27.06 3rd Qu.: 0.26
## Max. : 83.37 Max. : 180.00 Max. :3015.55
## NA's :200864
我so2df.Rda 在这里保存了这个数据帧(22M).
num_points = dim(lat_bounds)[1]
pixel_corners <- data.frame(lat_bounds=as.vector(lat_bounds), lon_bounds=as.vector(lon_bounds))
# create id column by replicating pixel's id for each of the 4 corner points
pixel_corners$id <- rep(so2df$id, each=num_points)
saveRDS(pixel_corners, file="pixel_corners.Rda")
summary(pixel_corners)
## lat_bounds lon_bounds id
## Min. :-89.96 Min. :-180.00 Length:5839200
## 1st Qu.:-62.29 1st Qu.:-163.30 Class :character
## Median :-19.86 Median :-150.46 Mode :character
## Mean :-13.87 Mean : -90.72
## 3rd Qu.: 31.26 3rd Qu.: -27.06
## Max. : 83.40 Max. : 180.00
正如预期的那样,纬度/经度边界数据帧是值数据帧的四倍(每个像素/值四个点).
我pixel_corners.Rda 在这里保存了这个数据帧(26M).
合并数据帧
然后我通过id合并两个数据帧:
start_time <- Sys.time()
so2df <- merge(pixel_corners, so2df, by=c("id"))
time_taken <- Sys.time() - start_time
print(paste(dim(so2df)[1], "rows merged in", time_taken, "seconds"))
## [1] "5839200 rows merged in 42.4763631820679 seconds"
正如您所看到的,这是一个CPU密集型过程.我想知道如果我一次使用15个文件(全球覆盖)会发生什么.
绘制数据
现在我已经将像素角连接到像素值,我可以轻松地绘制它们.通常,我对轨道的特定区域感兴趣,所以我创建了一个函数,在绘制输入数据帧之前对其进行子集化:
PlotRegion <- function(so2df, latlon, title) {
# Plot the given dataset over a geographic region.
#
# Args:
# df: The dataset, should include the no2tc, lat, lon columns
# latlon: A vector of four values identifying the botton-left and top-right corners
# c(latmin, latmax, lonmin, lonmax)
# title: The plot title
# subset the data frame first
df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])
subtitle = paste("#Pixel =", dim(df_sub)[1], "- Data min =",
formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2), "max =",
formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))
ggplot(df_sub) +
geom_polygon(aes(y=lat_bounds, x=lon_bounds, fill=so2tc, group=id), alpha=0.8) +
borders('world', xlim=range(df_sub$lon), ylim=range(df_sub$lat),
colour='gray20', size=.2) +
theme_light() +
theme(panel.ontop=TRUE, panel.background=element_blank()) +
scale_fill_distiller(palette='Spectral') +
coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
labs(title=title, subtitle=subtitle,
x="Longitude", y="Latitude",
fill=expression(DU))
}
然后我在感兴趣的区域调用我的函数,例如让我们看看在夏威夷发生的事情:
latlon = c(17.5, 22.5, -160, -154)
PlotRegion(so2df, latlon, expression(SO[2]~total~vertical~column))
它们是我的像素,看起来像是来自Mauna Loa的二氧化硫羽流.请暂时忽略负值.如您所见,像素区域朝着条带的边缘变化(不同的分箱方案).
我尝试使用ggmap在谷歌的地图上显示相同的情节:
PlotRegionMap <- function(so2df, latlon, title) {
# Plot the given dataset over a geographic region.
#
# Args:
# df: The dataset, should include the no2tc, lat, lon columns
# latlon: A vector of four values identifying the botton-left and top-right corners
# c(latmin, latmax, lonmin, lonmax)
# title: The plot title
# subset the data frame first
df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])
subtitle = paste("#Pixel =", dim(df_sub)[1], "Data min =", formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2),
"max =", formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))
base_map <- get_map(location = c(lon = (latlon[4]+latlon[3])/2, lat = (latlon[1]+latlon[2])/2), zoom = 7, maptype="terrain", color="bw")
ggmap(base_map, extent = "normal") +
geom_polygon(data=df_sub, aes(y=lat_bounds, x=lon_bounds,fill=so2tc, group=id), alpha=0.5) +
theme_light() +
theme(panel.ontop=TRUE, panel.background=element_blank()) +
scale_fill_distiller(palette='Spectral') +
coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
labs(title=title, subtitle=subtitle,
x="Longitude", y="Latitude",
fill=expression(DU))
}
这就是我得到的:
latlon = c(17.5, 22.5, -160, -154)
PlotRegionMap(so2df, latlon, expression(SO[2]~total~vertical~column))
问题
- 有没有更有效的方法来解决这个问题？我正在阅读关于
sf包的内容,我想知道我是否可以定义点的数据帧(值+中心像素坐标),并sf自动推断像素边界.这将使我不必依赖原始数据集中定义的纬度/经度边界,并且必须将它们与我的值合并.我可以接受朝向条带边缘的过渡区域的精度损失,否则网格非常规则,每个像素大小为3.5x7平方公里. - 是否可以通过聚合相邻像素来将我的数据重新网格化到常规网格(如何？),从而提高性能？我正在考虑使用这个
raster包,据我所知,它需要在常规网格上使用数据.这应该是有用的全球范围(例如欧洲的情节),我不需要绘制单个像素 - 实际上,我甚至都看不到它们. - 在绘制谷歌地图时,是否需要重新投影我的数据？
[奖励化妆品问题]
- 有没有更优雅的方法将我的数据框子集在由其四个角点标识的区域上？
- 如何更改色标以使较高的值相对于较低的值突出？我对日志规模经验不佳,效果不佳.
我想data.table这里可能会有帮助。合并几乎是立即的。
“1.24507117271423 秒内合并了 5839200 行”
library(data.table)
pixel_cornersDT <- as.data.table(pixel_corners)
so2dfDT <- as.data.table(so2df)
setkey(pixel_cornersDT, id)
setkey(so2dfDT, id)
so2dfDT <- merge(pixel_cornersDT, so2dfDT, by=c("id"), all = TRUE)
将数据放在 a 中data.table，绘图函数中的子集化也会快得多。
- 问题 1/2/4：
raster我认为使用or不会使该过程更快sf，但您可以尝试使用函数rasterFromXYZ()or st_make_grid()。但大部分时间将花费在转换为栅格/SF 对象上，因为您必须转换整个数据集。
我建议进行所有数据处理，包括data.table裁剪，然后您可以切换到光栅/SF 对象以进行绘图。
- 问题3：
谷歌图显示正确，但您指定了黑白地图，并且其覆盖有“光栅”，因此您不会看到很多。您可以将底图更改为卫星背景。
base_map <- get_map(location = c(lon = (latlon[4]+latlon[3])/2, lat = (latlon[1]+latlon[2])/2),
zoom = 7, maptype="satellite")
- 问题5：
您可以使用包rescale中的功能scales。我在下面提供了两个选项。第一个（未注释）将分位数作为中断，其他中断是单独定义的。我不会使用对数转换（trans- 参数），因为您会创建 NA 值，因为您也有负值。
ggplot(df_sub) +
geom_polygon(aes(y=lat_bounds, x=lon_bounds, fill=so2tc, group=id), alpha=0.8) +
borders('world', xlim=range(df_sub$lon), ylim=range(df_sub$lat),
colour='gray20', size=.2) +
theme_light() +
theme(panel.ontop=TRUE, panel.background=element_blank()) +
# scale_fill_distiller(palette='Spectral', type="seq", trans = "log2") +
scale_fill_distiller(palette = "Spectral",
# values = scales::rescale(quantile(df_sub$so2tc), c(0,1))) +
values = scales::rescale(c(-3,0,1,5), c(0,1))) +
coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
labs(title=title, subtitle=subtitle,
x="Longitude", y="Latitude",
fill=expression(DU))
整个过程现在对我来说大约需要8 秒，包括没有背景地图的绘图，尽管地图渲染也需要额外的 1-2 秒。