NCL批量绘制CMA台风路径（基于地球影像）

小其其格 · 发表于 2020-10-6 20:08:58

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由小其其格于 2020-10-7 15:31 编辑

---------------------------------- 我是分割线 ------------------------------------------------------------------------------
最新效果图：

---------------------------------- 我是分割线 ------------------------------------------------------------------------------

看到@夏朗的芒果大神的帖子：http://bbs.06climate.com/forum.php?mod=viewthread&tid=57969，画出来的台风路径效果太棒啦，所以也参照他的帖子，基于地球影像和中国气象局热带气旋资料中心提供的数据（http://tcdata.typhoon.org.cn/zjljsjj_zlhq.html），借助NCL批量绘制台风路径。

---------------------------------- 我是分割线 ---------------------------------------------------------------------------------
step1. 根据官网方法，将 EarthMap_2500x1250.jpg 转换成 EarthMap_2500x1250.nc
step2. 根据官网方法http://www.ncl.ucar.edu/Applications/Scripts/topo_9.ncl ，复制其 recreate_jpeg_image 函数，并进行相应修改
step3. 台风路径绘制脚本
---------------------------------- 我是分割线 ---------------------------------------------------------------------------------

子函数部分为：

;----------------------------------------------------------------------
; function recreate_jpeg_image
;----------------------------------------------------------------------
undef("recreate_jpeg_image")
function recreate_jpeg_image(wks,minlat,maxlat,minlon,maxlon,tcname_cma,tcname_cma2)
begin
jpg_filename = "EarthMap_2500x1250.jpg"
nc_filename = "EarthMap_2500x1250.nc"
;--You could use a system call to do the NetCDF conversion
; cmd = "gdal_translate -ot Int16 -of netCDF " + jpg_filename + \
; " " + nc_filename
; system(cmd)
;---Read the three bands of data
f = addfile(nc_filename,"r")
Band1 = where(f->Band1.gt.255, 255, f->Band1) ; red channel
Band2 = where(f->Band2.gt.255, 255, f->Band2) ; green channel
Band3 = where(f->Band3.gt.255, 255, f->Band3) ; blue channel
band_dims = dimsizes(Band3)
nlat = band_dims(0)
nlon = band_dims(1)
print("dimensions of image = " + nlat + " x " + nlon)
;
; Add lat/lon data so we can overlay on a map, and/or
; overlay contours. We know the image is global,
; cylindrical equidistant, and centered about lon=0.
;
lat = fspan( -90, 90,nlat)
lon = fspan(-180,180,nlon)
lat@units = "degrees_north"
lon@units = "degrees_east"
Band1!0 = "lat"
Band1!1 = "lon"
Band2!0 = "lat"
Band2!1 = "lon"
Band3!0 = "lat"
Band3!1 = "lon"
Band1&lat = lat
Band1&lon = lon
Band2&lat = lat
Band2&lon = lon
Band3&lat = lat
Band3&lon = lon
res = True
; res@vpWidthF = 0.85 ; Force image to fill screen.
; res@vpHeightF = 0.42
res@gsnMaximize = True
res@gsnFrame = False ; Don't draw or advance
res@gsnDraw = False ; frame yet.
res@cnFillOn = True
res@cnFillMode = "RasterFill" ; Raster fill can be faster
res@cnFillDrawOrder = "PreDraw"
; res@cnRasterCellSizeF = 0.0003
res@cnLevelSelectionMode = "EqualSpacedLevels"
res@cnMaxLevelCount = 254
res@cnFillBackgroundColor = (/ 1., 1., 1., 1./)
res@cnLinesOn = False ; Turn off contour lines .
res@cnLineLabelsOn = False ; Turn off contour labels
res@cnInfoLabelOn = False ; Turn off info label
res@lbLabelBarOn = False ; Turn off labelbar
res@gsnRightString = "" ; Turn off subtitles
;res@gsnLeftString = ""
;res@gsnLeftString = "track_"+tcname_cma+"_"+tcname_cma2 ; add the gsn/ti titles
; res@gsnTickMarksOn = False ; Turn off tickmarks
res@pmTickMarkDisplayMode = "Always"
;;---- add the title
res@tiMainString = "TRACK_"+tcname_cma+"_"+tcname_cma2 ; add titles
;res@tiMainString = "Main title left-justified"
;res@tiMainJust = "CenterLeft"
;res@tiMainPosition = "Left"
;---Construct RGBA colormaps...
ramp = fspan(0., 1., 255)
reds = new((/255, 4/), float)
greens = new((/255, 4/), float)
blues = new((/255, 4/), float)
reds = 0
greens = 0
blues = 0
reds(:,0) = ramp
greens(:,1) = ramp
blues(:,2) = ramp
; The red contour map is plotted fully opaque; the green and blue
; are plotted completely transparent. When overlain, the colors
; combine (rather magically).
reds(:,3) = 1.
greens(:,3) = 0
blues(:,3) = 0
res@cnFillColors = greens
greenMap = gsn_csm_contour(wks, Band2, res)
res@cnFillColors = blues
blueMap = gsn_csm_contour(wks, Band3, res)
;---This will be our base, so make it a map plot.
res@gsnAddCyclic = False
; res@mpShapeMode = "FreeAspect"
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "Earth..4"
res@mpOutlineSpecifiers = (/"China:Provinces","Taiwan"/)
res@mpNationalLineColor = "yellow"
res@mpProvincialLineColor = "yellow"
res@mpNationalLineThicknessF = 9.0
res@mpProvincialLineThicknessF = 5.0
res@mpGeophysicalLineColor = "white"
res@mpGeophysicalLineThicknessF = 1.5
res@mpOutlineDrawOrder = "Draw"
res@mpFillOn = False
res@mpGridAndLimbOn = True
res@mpGridLineThicknessF = 0.8
res@mpGridLineColor = "gray"
res@mpGridLatSpacingF = 5.
res@mpGridLonSpacingF = 5.
res@gsnMajorLatSpacing = 5.
res@gsnMajorLonSpacing = 5.
;---Zoom in on area of interest
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
res@cnFillColors = reds
redMap = gsn_csm_contour_map(wks, Band1, res)
;---Overlay everything to create the topo map
overlay(redMap, greenMap)
overlay(redMap, blueMap)
return(redMap)
end

复制代码

主函数部分，读取CMA提供的台风数据并批量绘制台风路径（增加了标题、台风标签以及图例，标注00时对应日期和放大对应圆圈）：

--------------------------------- 我是分割线 ----------------------------------------------------------------------------------

效果图：

---------------------------------- 我是分割线 ------------------------------------------------------------------------------

底图图像和测试数据附件下载：

游客，如果您要查看本帖隐藏内容请回复

freekiller · 发表于 2020-10-6 21:52:37

赞。谢谢分享

GYP123 · 发表于 2020-10-7 09:51:11

赞，谢谢分享

t_imba · 发表于 2020-10-8 06:18:24

谢谢大佬分享

小其其格 · 发表于 2020-10-8 07:19:19

t_imba 发表于 2020-10-8 06:18
谢谢大佬分享

共同进步

Zhefan · 发表于 2020-10-8 16:15:43

不错不错不错不错不错不错不错

Beckyeahyeah · 发表于 2020-10-10 08:31:52

1111111111111111111

花容天下重莲 · 发表于 2020-10-10 09:13:58

太厉害了，多向大佬学习。

鱼丸粗面 · 发表于 2020-10-10 11:42:20

谢谢大佬分享

cycle · 发表于 2020-10-11 09:24:13

谢谢大佬分享

		自动登录	找回密码
密码			立即注册

[作图] NCL批量绘制CMA台风路径（基于地球影像）

登录后查看更多精彩内容~