爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 8356|回复: 9

[作图] 关于add China map的问题

[复制链接]

新浪微博达人勋

发表于 2019-1-17 15:31:13 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
求各位大神指点
画panel图时使用add_china_map添加地图 只能加到最后一张图上 不知道是哪里的问题?

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
begin
;==========load types file============
f01                      = asciiread ("F:/cygwin/home/user01/天气型/ST.txt" , (/4681,1/),"float")
types                    = new((/4681,1/),float)
types@_FillValue         = 99999
types                    = f01
;=====================================
;==========load NCEP surface temperture file ==========
f02                      = addfile("J:/ncepdaily/surfterm19612013.nc", "r")
ST                       = f02->ST
;======================================================
;==========split weather types==========
ST01                     = new((/801,73,144/), float)
ST02                     = new((/653,73,144/), float)
ST03                     = new((/761,73,144/), float)
ST04                     = new((/611,73,144/), float)
ST05                     = new((/1044,73,144/), float)
ST06                     = new((/811,73,144/), float)
j=0
k=0
m=0
n=0
h=0
l=0
do i=0,4680
    if (types(i,0).eq.1) then
        ST01(j,:,:) = ST(i,:,:)
        j=j+1
    end if
    if (types(i,0).eq.2) then
        ST02(k,:,:) = ST(i,:,:)
        k=k+1
    end if
    if (types(i,0).eq.3) then
        ST03(m,:,:) = ST(i,:,:)
        m=m+1
    end if  
    if (types(i,0).eq.4) then
        ST04(n,:,:) = ST(i,:,:)
        n=n+1
    end if
    if (types(i,0).eq.5) then
        ST05(h,:,:) = ST(i,:,:)
        h=h+1
    end if  
    if (types(i,0).eq.6) then
        ST06(l,:,:) = ST(i,:,:)
        l=l+1
    end if
end do
ST001                     = dim_avg_n(ST01,0)
ST002                     = dim_avg_n(ST02,0)
ST003                     = dim_avg_n(ST03,0)
ST004                     = dim_avg_n(ST04,0)
ST005                     = dim_avg_n(ST05,0)
ST006                     = dim_avg_n(ST06,0)
ST001_re                  = new((/73,144/), float)
ST002_re                  = new((/73,144/), float)
ST003_re                  = new((/73,144/), float)
ST004_re                  = new((/73,144/), float)
ST005_re                  = new((/73,144/), float)
ST006_re                  = new((/73,144/), float)
do i = 0,72
    do j = 0,71
    ST001_re(i,j)=ST001(72-i,j+72)
    ST002_re(i,j)=ST002(72-i,j+72)
    ST003_re(i,j)=ST003(72-i,j+72)
    ST004_re(i,j)=ST004(72-i,j+72)
    ST005_re(i,j)=ST005(72-i,j+72)
    ST006_re(i,j)=ST006(72-i,j+72)
    end do
    do k = 72,143
    ST001_re(i,k)=ST001(72-i,k-72)
    ST002_re(i,k)=ST002(72-i,k-72)
    ST003_re(i,k)=ST003(72-i,k-72)
    ST004_re(i,k)=ST004(72-i,k-72)
    ST005_re(i,k)=ST005(72-i,k-72)
    ST006_re(i,k)=ST006(72-i,k-72)
    end do
end do
ST_RE = new((/6,73,144/),float)
ST_RE(0,:,:) = ST001_re
ST_RE(1,:,:) = ST002_re
ST_RE(2,:,:) = ST003_re
ST_RE(3,:,:) = ST004_re
ST_RE(4,:,:) = ST005_re
ST_RE(5,:,:) = ST006_re
;=====================================
;==========draw these pictures========
wks = gsn_open_wks("png","ST1")  
gsn_define_colormap(wks, "MPL_YlOrRd")
plot1 = new(6,graphic)
res                         = True            
res@tiMainString            = ""
res@gsnMaximize             = True
res@gsnDraw                 = False
res@gsnFrame                = False      
mpres                       = res
mpres@mpMinLatF               = 37.                        
mpres@mpMaxLatF               = 42.
mpres@mpMaxLonF               = 120
mpres@mpMinLonF               = 113
mpres@mpFillOn                = False
mpres@mpOutlineOn             = True  ; Use outlines from shapefile
mpres@mpDataBaseVersion       = "MediumRes"
mpres@mpDataSetName           = "F:/cygwin/app/ncl/lib/ncarg/nclscripts/data/database/Earth..4"
mpres@mpAreaMaskingOn         = True
mpres@mpMaskAreaSpecifiers    = (/"Beijing"/)
mpres@mpLandFillColor         = "white"
mpres@mpInlandWaterFillColor  = "white"
mpres@mpOceanFillColor        = "white"
mpres@mpOutlineBoundarySets   = (/"China","China:Provinces"/)
mpres@mpNationalLineColor        = "black"
mpres@mpProvincialLineColor      = "black"
mpres@mpGeophysicalLineColor     = "black"
mpres@mpNationalLineThicknessF   = 2
mpres@mpProvincialLineThicknessF = 1
mpres@mpGeophysicalLineThicknessF =3
mpres@tmXTOn                  = False
mpres@tmYROn                  = False
mpres@tmXBMode                = "Explicit"
mpres@tmXBValues              = (/114,116,118,120/)
mpres@tmXBLabels              = (/"114E","116E","118E","120E"/)
mpres@tmXBLabelFontHeightF    = 0.025
mpres@tmXBLabelFontThicknessF = 2  
mpres@tmYLLabelsOn            = True
mpres@tmYLMode                = "Explicit"
mpres@tmYLValues              = (/37,39,41/)
mpres@tmYLLabels              = (/"37N","39N","41N"/)
mpres@tmYLLabelFontHeightF    = 0.025
mpres@tmYLLabelFontThicknessF = 2  
mpres@tmXBMajorThicknessF     = 3
mpres@tmYLMajorThicknessF     = 3
mpres@tmXBMinorThicknessF     = 0
mpres@tmYLMinorThicknessF     = 0

mpres@cnFillOn                = True               
mpres@cnLinesOn               = True
mpres@cnLineLabelsOn          = False
mpres@cnLineDashSegLenF          = 0.18         
mpres@cnLineLabelInterval        = 1            
mpres@cnLineLabelPlacementMode   = "constant"
mpres@cnLineLabelFontHeightF     = 0.008
mpres@cnLevelSelectionMode    = "ManualLevels"
mpres@cnMinLevelValF          = 256
mpres@cnMaxLevelValF          = 282
mpres@cnLevelSpacingF         = 2
mpres@cnLineThicknessF        = 0
mpres@lbLabelsOn              = False
mpres@lbLabelBarOn            = False
mpres@tmXBLabelsOn =True

cnres                       = True
cnres@china                 = True       ;draw china map or not
cnres@river                 = True       ;draw changjiang&huanghe or not
cnres@province              = True       ;draw province boundary or not
cnres@nanhai                = False       ;draw nanhai or not
cnres@diqu                  = False       ; draw diqujie or not
do i=0,5
plot1(i) = gsn_csm_contour_map(wks,ST_RE(i,:,:),mpres)
chinamap = add_china_map(wks,plot1(i),cnres)
end do


pnlres = True ;
pnlres@gsnPanelDebug = True ;
pnlres@gsnPanelBottom = 0.05
pnlres@gsnPanelSave = True ; Save the state of the paneled plots so we can  
pnlres@gsnPanelLabelBar = True
pnlres@lbLabelFontHeightF = 0.015
pnlres@lbLabelFontThicknessF = 2
pnlres@lbBoxLineThicknessF= 2
pnlres@pmLabelBarWidthF   = 0.6
cnres                       = True
cnres@china                 = True       ;draw china map or not
cnres@river                 = True       ;draw changjiang&huanghe or not
cnres@province              = True       ;draw province boundary or not
cnres@nanhai                = False       ;draw nanhai or not
cnres@diqu                  = False       ; draw diqujie or not
gsn_panel(wks,plot1,(/3,2/),pnlres)  

draw(plot1)
frame(wks)
end
求大神指点!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-1-17 15:36:26 | 显示全部楼层
图片在这里
ST1.000001.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-1-17 15:54:29 | 显示全部楼层
chinamap = add_china_map(wks,plot1(i),cnres)
这里试一下chinamap也设成和plot1一样的数组,即
chinamap(i) = add_china_map(wks,plot1(i),cnres)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2019-1-17 16:19:31 | 显示全部楼层
croton 发表于 2019-1-17 15:54
chinamap = add_china_map(wks,plot1(i),cnres)
这里试一下chinamap也设成和plot1一样的数组,即
chinama ...

谢谢你的建议 使用chinamap(i)这种方法需要对chinamap进行定义 循环的结果会出错
我把这块拿出循环写成这样:
contour1 = add_china_map(wks,plot1(0),cnres)
contour2 = add_china_map(wks,plot1(1),cnres)
contour3 = add_china_map(wks,plot1(2),cnres)
contour4 = add_china_map(wks,plot1(3),cnres)
contour5 = add_china_map(wks,plot1(4),cnres)
contour6 = add_china_map(wks,plot1(5),cnres)
就可以达到先要的效果了。
ST1.000001.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-1-17 16:42:19 | 显示全部楼层
循环不会出错,只是楼主需要向plot1那样先定义一个图像数组,然后再循环。ncl不会自动设置数组,使用前都需要先定义。
chinamap = new(6,graphic)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-1-18 11:04:36 | 显示全部楼层
可以试一下楼上的方法 先定义chinamap = new(6,graphic) 然后在循环里写chinamap = add_china_map(wks,plot1(i),cnres) chinamap后不要加(i)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-11-22 15:26:20 | 显示全部楼层
不得不感慨一下万能的气象家园,遇到问题立马解决了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-11-24 11:25:13 | 显示全部楼层
因为也出现了楼主一样的问题 所以记录一下
前面楼里说的chinamap = new(6,graphic)再进行循环并不能解决问题,因为维度不匹配,输出看一下发现其维度是2172,也就是并非单个graphic,所以建立数组:
chinamap = new((/6,2172/),graphic)
再进行循环,成功
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-9-25 19:02:10 | 显示全部楼层
飛行船 发表于 2020-11-24 11:25
因为也出现了楼主一样的问题 所以记录一下
前面楼里说的chinamap = new(6,graphic)再进行循环并不能解决问 ...

请问再如何循环,谢谢您
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-9-25 19:28:09 | 显示全部楼层
甜酒果 发表于 2021-9-25 19:02
请问再如何循环,谢谢您

解决 :chinamap(i,:)        = add_china_map(wks,plot(i),cnres)  
感谢飛行船!!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表