问题描述
从 GSHHS 海岸线形状文件创建高分辨率土地掩膜
library(rgdal)
library(raster)
source("utils.R")
gshhs_dir = "C:/Users/mitch/Documents/R/Reef_Futures/gshhg-shp-2.3.7/GSHHS_shp"
合并除南极洲 (L1) 以外世界的全分辨率 ('f') 海岸线 南极洲 (L5) 的高分辨率 ('h') 海冰线 （由于在 (180,-90) 处缺少顶点，因此不对 L5 使用 'f'）
gshhs1 <- readOGR(file.path(gshhs_dir,"f"),"GSHHS_f_L1")
gshhs5 <- readOGR(file.path(gshhs_dir,"h"),"GSHHS_h_L5")
gshhs_full <- rbind(gshhs1,gshhs5)
将海岸线多边形栅格化为我们最终网格分辨率的 10 倍，1 = 陆地，0 = 海洋
new_rast <- raster(nrows = 43190,ncols = 86390,crs = CRS("+proj=longlat +datum=wgs84 +no_defs +ellps=wgs84 +towgs84=0,0"))
gshhs_rast <- rasterize(gshhs_full,new_rast,field = 1,background = 0,filename = "reeflandarea/gshhs_rast.Grd",datatype = "INT1U")
旋转到 (0,360) 经度范围
gshhs_rotate <- inv_rotate(gshhs_rast)
writeraster(gshhs_rotate,filename = "reeflandarea/gshhs_rotated.Grd",datatype = "INT1U")
从右边缘到左边缘'包裹'5个单元格，并在顶部和底部添加5个单元格 匹配最后一层的范围
extent_left <- extent(360 - 5*res(gshhs_rotate)[1],360,-90,90)
extent_right <- extent(0,360 - 5*res(gshhs_rotate)[1],90)
gshhs_left <- crop(gshhs_rotate,extent_left)
xmin(gshhs_left) <- xmin(gshhs_left) - 360
xmax(gshhs_left) <- 0
gshhs_right <- crop(gshhs_rotate,extent_right)
gshhs_top <- raster(nrows = 5,ncols = ncol(gshhs_rotate),vals = 0,ext = extent(xmin(gshhs_left),xmax(gshhs_right),90,90 + 5*res(gshhs_rotate)[2]))
gshhs_bottom <- raster(nrows = 5,vals = 1,-90 - 5*res(gshhs_rotate)[2],-90))
gshhs_wrap <- merge(gshhs_left,gshhs_right,gshhs_top,gshhs_bottom,filename = "reeflandarea/gshhs_wrap.Grd",datatype = "INT1U")
按因子 10 聚合以获得最终产品的分辨率 fun = min 以便单元格只有在所有基础值都是土地时才是土地
land_final <- aggregate(gshhs_wrap,fact=10,fun = min,datatype = "INT1U",filename = "reeflandarea/land_final.Grd")
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）