이수현
Reputation: 1
I'm trying to draw a weather chart using ncl, but I'm having a problem
I'm a student majoring in meteorology I'm trying to draw a weather chart using ncl, but I'm having two problems
The first question is that the High and Low labels are not drawn No matter how much I change the overlap setting, it won't be drawn on my picture
The second problem is that I use the Lambert conformal projection, so I used the tick marker code provided on the official website of ncl
But there's no space between tick marker and painting
Below is my code
1 load "add_map_tickmarks.ncl"
2
3 files := systemfunc("ls /home/shlee/2024_KIAPS/KIM_4.0.07/ne090np3_REAL_CASE/6hour_data/")
4
5 nfile = dimsizes(files)
6 dates = new(nfile, string)
7
8 explist = new(nfile, string)
9
10 do j = 0, nfile-1
11 dates(j) = str_get_field(files(j), 2, "-")
12 dates(j) = str_get_field(dates(j), 1, ".")
13 end do
14
15 ;------------------------------------------------------------------------------
16 ; resource
17 ;------------------------------------------------------------------------------
18
19 res = True
20
21 res@gsnMaximize = True
22
23 ;--- contour -------------------------------------------------------------------
24 ;res@cnHighLowLabelOverlapMode = "OmitOverVPAndHL"
25
26
27 res@cnLabelDrawOrder = "PostDraw"
28
29 res@cnFillOn = False
30 res@cnLinesOn = True
31
32 res@cnLevelSelectionMode = "ExplicitLevels"
33 res@cnLevels = ispan(4700, 6200, 30)
34
35 res@cnLineColor = "blue"
36 res@cnLineLabelFontColor = "blue"
37 res@cnLineThicknessF = 1.15
38 res@cnLineLabelInterval = 2
39 res@cnLineLabelPlacementMode = "constant"
40
41 res@cnInfoLabelOn = False
42
43 res@gsnFrame = False
44 res@gsnDraw = False
45
46 res@mpDataBaseVersion = "MediumRes"
47 res@mpDataSetName = "Earth..1"
48 ;--- string -------------------------------------------------------------------
49 res@gsnRightString = ""
50 res@gsnLeftString = ""
51
52 ;--- map plot -----------------------------------------------------------------
53 res@mpOutlineOn = True
54 res@mpOutlineBoundarySets = "AllBoundaries"
55
56 res@mpLimitMode = "LatLon"
57
58 res@mpMinLonF = 77.0
59 res@mpMaxLonF = 160.0
60 res@mpMinLatF = 12.5
61 res@mpMaxLatF = 60.0
62
63 res@mpOceanFillColor = "lightcyan"
64 res@mpLandFillColor = "lightyellow1"
65 res@mpGeophysicalLineThicknessF = 2
66
67 res@mpProjection = "LambertConformal"
68 res@mpLambertParallel1F = 35.
69 res@mpLambertParallel2F = 45.
70 res@mpLambertMeridianF = 125.
71
72
73 res@mpGridAndLimbOn = True
74 res@mpGridLineThicknessF = 0.8
75 res@mpGridSpacingF = 10.0
76
77 ;--- figure ratio -------------------------------------------------------------
78 res@mpShapeMode = "FreeAspect"
79 res@vpWidthF = 0.9
80 res@vpHeightF = 0.65
81
82 res@pmTickMarkDisplayMode = "Always"
83
84
85 ;--- Temperature contour resource --------------------------------------------
86 temres = True
87
88 temres@gsnMaximize = True
89 temres@gsnDraw = False
90 temres@gsnFrame = False
91
92 temres@cnFillOn = True
93 temres@cnLinesOn = False
94 temres@cnFillPallete = "MPL_rainbow"
95
96 ;temres@cnLevelSelectionMode = "ExplicitLevels"
97 ;temres@cnLevels = ispan(4700, 6200, 30)
98
99 res@cnInfoLabelOn = False
100
101 ;------------------------------------------------------------------------------
102 ; plot
103 ;------------------------------------------------------------------------------
104
105 ;do i = 0, nfile-1
106 do i = 0, 40
107 wks = gsn_open_wks("png", "./REAL_FCST_hgt_rain_" + dates(i))
108
109 dash_pattern = (/"$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_$_"/)
110 new_dash = NhlNewDashPattern(wks, dash_pattern)
111
112 res@mpGridLineDashPattern = new_dash
113
114 res@cnHighLabelsOn = True
115 res@cnHighLabelFontHeightF = 0.03
116 res@cnHighLabelBackgroundColor = -1
117 res@cnHighLabelFontColor = "blue"
118
119 res@cnLowLabelsOn = True
120 res@cnLowLabelFontHeightF = 0.015
121 res@cnLowLabelBackgroundColor = -1
122 res@cnLowLabelFontColor = "red"
123 ncfile = addfile("/home/shlee/2024_KIAPS/KIM_4.0.07/ne090np3_REAL_CASE/post/REAL_post_" \
124 + dates(i) + ".nc", "r")
125
126 hgt = ncfile->hgt(0,13,:,:)
127 t = ncfile->T(0,13,:,:)
128
129 plot = gsn_csm_contour_map(wks, hgt, res)
130
131 ;--- tick marker resource -----------------------------------------------------
132 tmres = True
133
134 tmres@tmXBValues = ispan(90, 150, 10)
135 tmres@tmXTValues = (/ 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, -170, -160/)
136 tmres@tmYLValues = ispan(10, 40, 10)
137 tmres@tmYRValues = ispan(10, 50, 10)
138
139 tmYUseLeft = False
140 tmXUseBottom = False
141
142 tmres@tmXBLabelFontHeightF = 0.008
143 res@tmXBMajorLengthF = 0.001
144 res@tmXTMajorLengthF = 0.001
145 res@tmYLMajorLengthF = 0.001
146 res@tmYRMajorLengthF = 0.001
147
148 tmres@tmXBMajorLengthF = 0.001
149 tmres@tmXTMajorLengthF = 0.001
150 tmres@tmYLMajorLengthF = 0.001
151 tmres@tmYRMajorLengthF = 0.001
152 map = add_map_tickmarks2(wks, plot, tmres)
153
154 maximize_output(wks, True)
155
156
157 end do
1 undef("add_map_tickmarks2")
2 function add_map_tickmarks2(wks,plot,res)
3 local res2, bres, vpx, vpy, vpw, vph, xndc, yndc, npts, n, j, nlat, \
4 nlon, delta, bot_lon, top_lon, lft_lat, rgt_lat, xblabels, xbvalues, \
5 xtlabels, xtvalues, yllabels, ylvalues, yrlabels, yrvalues, xfix, \
6 xlat, xlon, yfix, annoid, anno_str
7 begin
8 ;---Make a copy of the original resource list.
9 res2 = res
10
11 ;---Retrieve edges of plot in NDC space.
12 getvalues plot
13 "vpXF" : vpx
14 "vpYF" : vpy
15 "vpWidthF" : vpw
16 "vpHeightF" : vph
17 end getvalues
18
19 ;---Turn off tickmarks associated with map. We want to add our own.
20 setvalues plot
21 "pmTickMarkDisplayMode" : "Never"
22 end setvalues
23
24 ;---Initialize resources for tickmark plot. User shouldn't change these.
25 bres = True
26 bres@vpXF = vpx
27 bres@vpYF = vpy
28 bres@vpWidthF = vpw
29 bres@vpHeightF = vph
30 bres@trXMinF = vpx
31 bres@trXMaxF = vpx + vpw
32 bres@trYMinF = vpy - vph
33 bres@trYMaxF = vpy
34 bres@tmEqualizeXYSizes = True
35
36 ;---This resource the user can change in main code if desired.
37 bres@gsnTickMarksPointOutward = get_res_value(res2,"gsnTickMarksPointOutward",True)
38
39 ;
40 ; NDC Points to scan on X and Y axes. These arrays will be used to
41 ; find the closest NDC pair that gets us close to the location where
42 ; we want a tickmark.
43 ;
44 npts = 100000 ; Increase to get closer match for tickmarks
45 xndc = fspan(vpx,vpx+vpw,npts)
46 yndc = fspan(vpy-vph,vpy,npts)
47
48 n = dimsizes(yndc)
49 xfix = new(n,float)
50 yfix = new(n,float)
51 xlon = new(n,float)
52 xlat = new(n,float)
53 delta = 0.001
54
55 ;---Left axis tickmarks
56 if(isatt(res2,"tmYLValues")) then
57 lft_lat = get_res_value(res2,"tmYLValues",-1)
58 nlat = dimsizes(lft_lat)
59 ylvalues = new(nlat,float)
60 yllabels = new(nlat,string)
61 xfix = vpx + 0.0001 ; Just a smidge into the plot to make sure we don't
62 ; get missing values returned.
63 ;
64 ; Loop across each left latitude value that we want a tickmark for,
65 ; and try to find the closest X,Y NDC coordinate pair along this axis.
66 ;
67 NhlNDCToData(plot,xfix,yndc,xlon,xlat)
68 do j=0,dimsizes(lft_lat)-1
69 NhlNDCToData(plot,xfix,yndc,xlon,xlat-0.001)
70 ii = minind(fabs(xlat-lft_lat(j)))
71 if(.not.any(ismissing(ii)).and.fabs(xlat(ii)-lft_lat(j)).le.delta)
72 yllabels(j) = fabs(lft_lat(j)) + ""
73 ylvalues(j) = yndc(ii(0))
74 if(lft_lat(j).lt.0) then
75 yllabels(j) = yllabels(j) + "~S~o~N~S"
76 end if
77 if(lft_lat(j).gt.0) then
78 yllabels(j) = yllabels(j) + "~S~o~N~N"
79 end if
80 end if
81 delete(ii)
82 end do
83 bres@tmYLMode = "Explicit"
84 bres@tmYLValues = ylvalues
85 bres@tmYLLabels = get_res_value(res2,"tmYLLabels",yllabels)
86 else
87 bres@tmYLOn = False
88 bres@tmYLLabelsOn = False
89 end if
90
91 ;---Right axis tickmarks
92 if(isatt(res2,"tmYRValues")) then
93 rgt_lat = get_res_value(res2,"tmYRValues",-1)
94 nlat = dimsizes(rgt_lat)
95 yrvalues = new(nlat,float)
96 yrlabels = new(nlat,string)
97
98 xfix = vpx + vpw - 0.0001 ; Just a smidge into the plot to make sure we don't
99 ; get missing values returned.
100 ;
101 ; Loop across each right latitude value that we want a tickmark for,
102 ; and try to find the closest X,Y NDC coordinate pair along this axis.
103 ;
104 do j=0,dimsizes(rgt_lat)-1
105 NhlNDCToData(plot,xfix,yndc,xlon,xlat)
106 ii = minind(fabs(xlat-rgt_lat(j)))
107 if(.not.any(ismissing(ii)).and.fabs(xlat(ii)-rgt_lat(j)).le.delta)
108 yrlabels(j) = fabs(rgt_lat(j)) + ""
109 yrvalues(j) = yndc(ii(0))
110 if(rgt_lat(j).lt.0) then
111 yrlabels(j) = yrlabels(j) + "~S~o~N~S"
112 end if
113 if(rgt_lat(j).gt.0) then
114 yrlabels(j) = yrlabels(j) + "~S~o~N~N"
115 end if
116 end if
117 delete(ii)
118 end do
119 bres@tmYROn = True
120 bres@tmYRLabelsOn = True
121 bres@tmYUseLeft = False
122 bres@tmYRMode = "Explicit"
123 bres@tmYRValues = yrvalues
124 bres@tmYRLabels = get_res_value(res2,"tmYRLabels",yrlabels)
125 else
126 bres@tmYUseLeft = False
127 bres@tmYROn = False
128 bres@tmYRLabelsOn = False
129 end if
130
131 ;---Top axis tickmarks
132 if(isatt(res2,"tmXTValues")) then
133 top_lon = get_res_value(res2,"tmXTValues",-1)
134 nlon = dimsizes(top_lon)
135 xtvalues = new(nlon,float)
136 xtlabels = new(nlon,string)
137
138 yfix = vpy - 0.0001 ; Just a smidge into the plot to make sure we don't
139 ; get missing values returned.
140 ;
141 ; Loop across each top longitude value that we want a tickmark for,
142 ; and try to find the closest X,Y NDC coordinate pair along this axis.
143 ;
144 do j=0,dimsizes(top_lon)-1
145 NhlNDCToData(plot,xndc,yfix,xlon,xlat)
146 ii = minind(fabs(xlon-top_lon(j)))
147 if(.not.any(ismissing(ii)).and.fabs(xlon(ii)-top_lon(j)).le.delta)
148 xtlabels(j) = fabs(top_lon(j)) + ""
149 xtvalues(j) = xndc(ii(0))
150 if(top_lon(j).lt.0) then
151 xtlabels(j) = xtlabels(j) + "~S~o~N~W"
152 end if
153 if(top_lon(j).gt.0) then
154 xtlabels(j) = xtlabels(j) + "~S~o~N~E"
155 end if
156 end if
157 delete(ii)
158 end do
159 bres@tmXTOn = True
160 bres@tmXTLabelsOn = True
161 bres@tmXUseBottom = False
162 bres@tmXTMode = "Explicit"
163 bres@tmXTValues = xtvalues
164 bres@tmXTLabels = get_res_value(res2,"tmXTLabels",xtlabels)
165 else
166 bres@tmXUseBottom = False
167 bres@tmXTOn = False
168 bres@tmXTLabelsOn = False
169 end if
170
171 ;---Bottom axis tickmarks
172 if(isatt(res2,"tmXBValues")) then
173 bot_lon = get_res_value(res2,"tmXBValues",-1)
174 nlon = dimsizes(bot_lon)
175 xbvalues = new(nlon,float)
176 xblabels = new(nlon,string)
177
178 yfix = vpy-vph + 0.0001 ; Just a smidge into the plot to make sure
179 ; we don't get missing values returned.
180 ;
181 ; Loop across each bottom longitude value that we want a tickmark for,
182 ; and try to find the closest X,Y NDC coordinate pair along this axis.
183 ;
184 do j=0,dimsizes(bot_lon)-1
185 NhlNDCToData(plot,xndc,yfix,xlon,xlat)
186 ii = minind(fabs(xlon-bot_lon(j)))
187 if(.not.any(ismissing(ii)).and.fabs(xlon(ii)-bot_lon(j)).le.delta)
188 xblabels(j) = fabs(bot_lon(j)) + ""
189 xbvalues(j) = xndc(ii(0))
190 if(bot_lon(j).lt.0) then
191 xblabels(j) = xblabels(j) + "~S~o~N~W"
192 end if
193 if(bot_lon(j).gt.0) then
194 xblabels(j) = xblabels(j) + "~S~o~N~E"
195 end if
196 end if
197 delete(ii)
198 end do
199 bres@tmXBMode = "Explicit"
200 bres@tmXBValues = xbvalues
201 bres@tmXBLabels = get_res_value(res2,"tmXBLabels",xblabels)
202 else
203 bres@tmXBOn = False
204 bres@tmXBLabelsOn = False
205 end if
206
207 ;
208 ; Now that we are done figuring out where to put tickmarks, and
209 ; what labels to use, get any "tm" resources that might have been
210 ; set by the user, and create a blank plot with thes new tickmarks.
211 ;
212
213 ;---Get rest of user resources that were set with "tm".
214 bres = get_res_eq(res2,"tm")
215
216 bres = True ; Above call will set bres to True if no "tm" resources, so
217 ; make sure it is True still.
218
219 bres@gsnDraw = False
220 bres@gsnFrame = False
221
222 ;
223 ; Create blank plot with new tickmarks (don't use gsn_csm_blank_plot,
224 ; because it wants to scale the size of your X and Y axes.)
225 ;
226 blank = gsn_blank_plot(wks,bres)
227
228 ;
229 ; Attach new tickmarks to original plot. This will allow resizing
230 ; if desired. The default is to attach one plot to the center of
231 ; the other one. These two plots are already the same size.
232 ;
233 annoid = gsn_add_annotation(plot,blank,False)
234
235 ;
236 ; Be sure to return the annotation id, otherwise the
237 ; tickmarks will disappear.
238 ;
239 anno_str = unique_string("annoid")
240 plot@$anno_str$ = annoid
241
242 return(plot)
243 end
Upvotes: 0
Views: 22