基于地球影像的台风路径绘制脚本，超级酷炫
本帖最后由 夏朗的芒果 于 2018-11-27 08:19 编辑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 的下载地址为
**** Hidden Message *****
子函数部分为
undef("recreate_jpeg_image")
function recreate_jpeg_image(wks,minlat,maxlat,minlon,maxlon)
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@gsnTickMarksOn= False ; Turn off tickmarks
res@pmTickMarkDisplayMode = "Always"
;---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
主脚本为
**** Hidden Message *****
台风路径见附件
这个看看学习下 楼主的想法很有创意，厉害了！ 楼主辛苦了。看看。 楼主厉害！学习了！ 这么厉害啊 厉害了~！！！来学习! 楼主厉害了！这么一画确实好看多了 楼主厉害了。来学习 紧跟楼主学python，不对，学ncl