基于地球影像的台风路径绘制脚本，超级酷炫

夏朗的芒果

2017.12.21 更新
替代了更加精细美观的底图（下载地址 https://visibleearth.nasa.gov/view_cat.php?categoryID=1484），同时加强了路径绘制效果，形成了终极魔改版



-----------------------------------------------------------------------


以下为原帖：
废话不多说，先上效果图，这是我小试牛刀绘制的今年1705号奥鹿路径。当然作用肯定不止于此，这里权作抛砖引玉。


引子：一次偶然的机遇，我在ncl官网上看到ncl可以将在jpg图片上再次绘图，具体见官网http://www.ncl.ucar.edu/Applications/topo.shtml 第9个案例。由此想到，如果在地球影像上叠加台风路径，岂不十分美妙？说干就干。

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

EarthMap_2500x1250.nc 的下载地址为

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

子函数部分为

1. undef("recreate_jpeg_image")
   
2. function recreate_jpeg_image(wks,minlat,maxlat,minlon,maxlon)
   
3. begin
   
4.   jpg_filename = "EarthMap_2500x1250.jpg"
   
5.   nc_filename       = "EarthMap_2500x1250.nc"
   
6. 
   
7. ;--You could use a system call to do the NetCDF conversion
   
8. ; cmd = "gdal_translate -ot Int16 -of netCDF " + jpg_filename + \
   
9. ;        " " + nc_filename
   
10. ; system(cmd)
    
11. 
    
12. ;---Read the three bands of data
    
13.   f     = addfile(nc_filename,"r")
    
14.   Band1 = where(f->Band1.gt.255, 255, f->Band1)  ; red channel
    
15.   Band2 = where(f->Band2.gt.255, 255, f->Band2)  ; green channel
    
16.   Band3 = where(f->Band3.gt.255, 255, f->Band3)  ; blue channel
    
17. 
    
18.   band_dims = dimsizes(Band3)
    
19.   nlat      = band_dims(0)
    
20.   nlon      = band_dims(1)
    
21.   print("dimensions of image = " + nlat + " x " + nlon)
    
22. 
    
23. ;
    
24. ; Add lat/lon data so we can overlay on a map, and/or
    
25. ; overlay contours. We know the image is global,
    
26. ; cylindrical equidistant, and centered about lon=0.
    
27. ;
    
28.   lat       = fspan( -90, 90,nlat)
    
29.   lon       = fspan(-180,180,nlon)
    
30.   lat@units = "degrees_north"
    
31.   lon@units = "degrees_east"
    
32. 
    
33.   Band1!0   = "lat"
    
34.   Band1!1   = "lon"
    
35.   Band2!0   = "lat"
    
36.   Band2!1   = "lon"
    
37.   Band3!0   = "lat"
    
38.   Band3!1   = "lon"
    
39.   Band1&lat = lat
    
40.   Band1&lon = lon
    
41.   Band2&lat = lat
    
42.   Band2&lon = lon
    
43.   Band3&lat = lat
    
44.   Band3&lon = lon
    
45. 
    
46.   res                 = True
    
47.   ; res@vpWidthF        = 0.85          ; Force image to fill screen.
    
48.   ; res@vpHeightF       = 0.42
    
49.   res@gsnMaximize     = True
    
50. 
    
51.   res@gsnFrame        = False        ; Don't draw or advance
    
52.   res@gsnDraw         = False        ; frame yet.
    
53. 
    
54.   res@cnFillOn        = True
    
55.   res@cnFillMode      = "RasterFill" ; Raster fill can be faster
    
56.   res@cnFillDrawOrder = "PreDraw"
    
57.   ; res@cnRasterCellSizeF = 0.0003
    
58. 
    
59.   res@cnLevelSelectionMode  = "EqualSpacedLevels"
    
60.   res@cnMaxLevelCount       = 254
    
61.   res@cnFillBackgroundColor = (/ 1., 1., 1., 1./)
    
62. 
    
63.   res@cnLinesOn       = False              ; Turn off contour lines      .
    
64.   res@cnLineLabelsOn  = False              ; Turn off contour labels
    
65.   res@cnInfoLabelOn   = False              ; Turn off info label
    
66.   res@lbLabelBarOn    = False              ; Turn off labelbar
    
67.   res@gsnRightString  = ""                 ; Turn off subtitles
    
68.   res@gsnLeftString   = ""
    
69.   ; res@gsnTickMarksOn  = False              ; Turn off tickmarks
    
70.   res@pmTickMarkDisplayMode    = "Always"
    
71. 
    
72. ;---Construct RGBA colormaps...
    
73.   ramp   = fspan(0., 1., 255)
    
74.   reds   = new((/255, 4/), float)
    
75.   greens = new((/255, 4/), float)
    
76.   blues  = new((/255, 4/), float)
    
77. 
    
78.   reds   = 0
    
79.   greens = 0
    
80.   blues  = 0
    
81. 
    
82.   reds(:,0)   = ramp
    
83.   greens(:,1) = ramp
    
84.   blues(:,2)  = ramp
    
85. 
    
86.   ; The red contour map is plotted fully opaque; the green and blue
    
87.   ; are plotted completely transparent. When overlain, the colors
    
88.   ; combine (rather magically).
    
89.   reds(:,3)   = 1.
    
90.   greens(:,3) = 0
    
91.   blues(:,3)  = 0
    
92. 
    
93.   res@cnFillColors = greens
    
94.   greenMap = gsn_csm_contour(wks, Band2, res)
    
95. 
    
96.   res@cnFillColors = blues
    
97.   blueMap = gsn_csm_contour(wks, Band3, res)
    
98. 
    
99. ;---This will be our base, so make it a map plot.
    
100.   res@gsnAddCyclic             = False
     
101.   ; res@mpShapeMode = "FreeAspect"
     
102.   res@mpDataBaseVersion        = "MediumRes"
     
103.   res@mpDataSetName              = "Earth..4"
     
104.   res@mpOutlineSpecifiers        = (/"China:Provinces","Taiwan"/)
     
105.   res@mpNationalLineColor = "yellow"
     
106.   res@mpProvincialLineColor = "yellow"
     
107.   res@mpNationalLineThicknessF = 9.0
     
108.   res@mpProvincialLineThicknessF = 5.0
     
109.   res@mpGeophysicalLineColor = "white"
     
110.   res@mpGeophysicalLineThicknessF = 1.5
     
111.   res@mpOutlineDrawOrder = "Draw"
     
112.   res@mpFillOn                 = False
     
113. 
     
114.   res@mpGridAndLimbOn   = True
     
115.   res@mpGridLineThicknessF = 0.8
     
116.   res@mpGridLineColor   = "gray"
     
117.   res@mpGridLatSpacingF     =  5.
     
118.   res@mpGridLonSpacingF     =  5.
     
119.   res@gsnMajorLatSpacing    =  5.
     
120.   res@gsnMajorLonSpacing    =  5.
     
121. 
     
122. ;---Zoom in on area of interest
     
123.   res@mpMinLatF                = minlat
     
124.   res@mpMaxLatF                = maxlat
     
125.   res@mpMinLonF                = minlon
     
126.   res@mpMaxLonF                = maxlon
     
127. 
     
128.   res@cnFillColors             = reds
     
129.   redMap = gsn_csm_contour_map(wks, Band1, res)
     
130. 
     
131. ;---Overlay everything to create the topo map
     
132.   overlay(redMap, greenMap)
     
133.   overlay(redMap, blueMap)
     
134. 
     
135.   return(redMap)
     
136. end

复制代码

主脚本为

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

台风路径见附件
201705.txt (3.49 KB, 下载次数: 143)

2017-12-20 00:03 上传

点击文件名下载附件

随缘

这个看看学习下

Soaring

楼主的想法很有创意，厉害了！

熙熙攘攘

楼主辛苦了。看看。

liuxiaoyue

楼主厉害！学习了！

rabin_xu

这么厉害啊

Quincy

厉害了~！！！来学习!

煮茶的水

楼主厉害了！这么一画确实好看多了

打破砂锅纹到底

楼主厉害了。来学习

打破砂锅纹到底

紧跟楼主学python，不对，学ncl