getmapbounds.f
Go to the documentation of this file.
1 c> \ingroup core
2 c> Subroutine to collect up the MAP boundaries
3 c> for use in creating NADCON 5
4 c>
5 c> This CAN BE different than the GRID boundaries
6 c> as such:
7 c>
8 c> GRID boundaries will be just four values (n/s/w/e) for any region
9 c>
10 c> MAP boundaries will allow multiple maps to be made and may or may
11 c> not align with the GRID boundaries. Used to allow for more
12 c> "close up" maps and such, without the need to screw up the
13 c> MAP boundaries.
14 c>
15 c> \param[in] mapflag Map Generation Flag
16 c> \param[in] maxplots
17 c> \param[in] region region to get map bounds
18 c> \param[out] nplots number of plots generated
19 c> \param[in] olddtm source datum
20 c> \param[in] newdtm target datum
21 c> \param[out] bw western bound of plot (Array of length maxplots)
22 c> \param[out] be eastern bound of plot (Array of length maxplots)
23 c> \param[out] bs southern bound of plot (Array of length maxplots)
24 c> \param[out] bn northern bound of plot (Array of length maxplots)
25 c> \param[out] jm (Array of length maxplots)
26 c> \param[out] b1 (Array of length maxplots)
27 c> \param[out] b2 (Array of length maxplots)
28 c> \param[out] fn (Array of length maxplots)
29 c> \param[out] lrv (Array of length maxplots)
30 c> \param[out] rv0x (Array of length maxplots)
31 c> \param[out] rv0y (Array of length maxplots)
32 c> \param[out] rl0y (Array of length maxplots)
33 c>
34 c> Version for NADCON 5
35 c> Built upon the original version used in GEOCON v2.0
36 c> Do not use with GEOCON v2.0
37 c>
38 c> Broken down into sub-subroutines to make it easier
39 c> to swap out when I make different choices.
40 c>
41 c>
42 c> ## Changelog
43 c>
44 c>
45 c> ### 2016 08 29:
46 c> Taking in olddtm and newdtm now, and adding
47 c> code to use that to filter out "Saint" regions in Alaska
48 c> when plotting transformations not supported in those regions.
49 c>
50 c> ### 2016 08 26:
51 c> Used actual mercator projection math to compute
52 c> the exact reference vector and label locations 1/2 inch and
53 c> 3/4 inch respectively below the S/W corner of the plot.
54 c>
55 c> ### 2016 07 21:
56 c> Two new columns added to "map.parameters", which
57 c> have the location of the reference vector. Return a logical "lrv"
58 c> as true if there is an optional special location for the reverence vector.
59 c> Return as false if not. If true, return lon/lat coords of ref vector
60 c> origin in rv0x/rv0y. If false, return zeros in those fields.
61 c>
62 c> Also, compute "jm" on the fly, ignoring what is in the table. All plots
63 c> will now be forced PORTRAIT and forced no wider than 6" and no taller
64 c> than 8", while maintaining proper X/Y ratios in a Mercator projection.
65 c> That means, make the biggest plot possible, with the right ratio, that
66 c> is neither wider than 6" nor taller than 8" and then, whatever the width
67 c> of that largest plot is -- return that width in the "jm" field.
68 c>
69
70  subroutine getmapbounds(mapflag,maxplots,region,nplots,
71  *olddtm,newdtm,
72  *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y)
74 c - 2016 08 29: Taking in olddtm and newdtm now, and adding
75 c - code to use that to filter out "Saint" regions in Alaska
76 c - when plotting transformations not supported in those regions.
77
78 c - 2016 08 26: Used actual mercator projection math to compute
79 c - the exact reference vector and label locations 1/2 inch and
80 c - 3/4 inch respectively below the S/W corner of the plot.
81
82 c - Updated 2016 07 21: Two new columns added to "map.parameters", which
83 c - have the location of the reference vector. Return a logical "lrv"
84 c - as true if there is an optional special location for the reverence vector.
85 c - Return as false if not. If true, return lon/lat coords of ref vector
86 c - origin in rv0x/rv0y. If false, return zeros in those fields.
87 c -
88 c - Also, compute "jm" on the fly, ignoring what is in the table. All plots
89 c - will now be forced PORTRAIT and forced no wider than 6" and no taller
90 c - than 8", while maintaining proper X/Y ratios in a Mercator projection.
91 c - That means, make the biggest plot possible, with the right ratio, that
92 c - is neither wider than 6" nor taller than 8" and then, whatever the width
93 c - of that largest plot is -- return that width in the "jm" field.
94
95 c
96 c - Subroutine to collect up the MAP boundaries
97 c - for use in creating NADCON 5
98 c
99 c - This CAN BE different than the GRID boundaries
100 c - as such:
101 c
102 c - GRID boundaries will be just four values (n/s/w/e) for any region
103 c
104 c - MAP boundaries will allow multiple maps to be made and may or may
105 c - not align with the GRID boundaries. Used to allow for more
106 c - "close up" maps and such, without the need to screw up the
107 c - MAP boundaries.
108
109 c - Input:
110 c region, maxplots, mapflag
111 c olddtm, newdtm
112 c - Output:
113 c number of plots, their bounds, their GMT scales and
114 c label deltas
115
116 c - Version for NADCON 5
117 c - Built upon the original version used in GEOCON v2.0
118 c - Do not use with GEOCON v2.0
119
120 c - Broken down into sub-subroutines to make it easier
121 c - to swap out when I make different choices.
122
123  character*10 region,olddtm,newdtm
124  real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots)
125  real*4 jm(maxplots)
126  real*4 b1(maxplots),b2(maxplots)
127  character*10 fn(maxplots)
128  character*1 mapflag
129
130
131 c - 2016 07 21:
132  logical lrv(maxplots)
133  real*8 rv0x(maxplots),rv0y(maxplots),rl0y(maxplots)
134  character*3 c3x,c3y
135  character*2 c2x,c2y
136  real*8 pi,d2r
137
138  character*200 card
139
140
141  pi = 2.d0*dasin(1.d0)
142  d2r = pi/180.d0
143
144  ifile = 3
145
146  open(ifile,
147  *file='Data/map.parameters',
148  *status='old',form='formatted')
149
150 c - Get numerical version of mapflag
152
153 c - Get the header out of the way
156
157 c - Loop over all map parameters, finding the ones relevant to
158 c - our region and mapflag.
159  i = 0
161  if(card( 1:10).ne.region)goto 1
163  if(iflag0.gt.iflag)goto 1
164
165 c - 2016 08 29:
166 c - Specific to the case where we are looking at region=alaska,
167 c - and NOT region=<some saint island>, then for the specific
168 c - case of transforming from NAD 27 forward, let's skip
169 c - the <saint islands>
171  * iflag.eq.2 .and.
172  * iflag0.eq.2 .and.
173  * (card(16:25).eq.'stpaul' .or.
174  * card(16:25).eq.'stgeorge' .or.
175  * card(16:25).eq.'stlawrence' .or.
176  * card(16:25).eq.'stmatthew') .and.
178
179  i = i + 1
180  fn(i) = trim(card( 16: 25))
194 c - 2016 07 21:
199
200  bw(i) = dble(iwd) + dble(iwm)/60.d0
201  be(i) = dble(ied) + dble(iem)/60.d0
202  if(ind.lt.0)inm = -inm
203  if(isd.lt.0)ism = -ism
204  bn(i) = dble(ind) + dble(inm)/60.d0
205  bs(i) = dble(isd) + dble(ism)/60.d0
206  jm(i) = xjm
207  b1(i) = dble(ib1d) + dble(ib1m)/60.d0
208  b2(i) = dble(ib2d) + dble(ib2m)/60.d0
209
210 c - 2016 07 21 (Optional reference vector locations):
211  if(c3x.eq.'---')then
212  lrv(i) = .false.
213  rv0x(i) = 0
214  rv0y(i) = 0
215  else
216  lrv(i) = .true.
221  rv0x(i) = dble(rv0xd) + dble(rv0xm)/60.d0
222  rv0y(i) = dble(rv0yd) + dble(rv0ym)/60.d0
223  endif
224
225 c - 2016 07 21 (Forcing the map to keep its proper Mercator-projected X/Y ratio,
226 c - while filling the maximum space that does not exceed 6" wide nor 8" high)
227  dx = (be(i) - bw(i)) * d2r
228
229  q1 = dtan(bn(i)*d2r)
230  q2 = 1.d0 / dcos(bn(i)*d2r) ! Secant
231  yn = log(q1 + q2)
232
233  q1 = dtan(bs(i)*d2r)
234  q2 = 1.d0 / dcos(bs(i)*d2r) ! Secant
235  ys = log(q1 + q2)
236
237  dy = yn - ys
238
239  ratio = dx/dy
240
241 c - Max height and width
242 c xmxht = 8
243  xmxht = 6
244  xmxwd = 6
245
246  if(xmxht*ratio .gt. xmxwd)then
247  jm(i) = xmxwd
248 c - And height = xmxwd / ratio...
249  else
250  jm(i) = xmxht*ratio
251 c - And height = xmxht
252  endif
253
254 c - 2016 08 26
255 c See DRU-12, p. 56-57
256 c - dy/dphi is linear for latitudes south of 65. All of our plots
257 c - have at most a 65 degree south border. So I can use that for
258 c - a pretty good approximation to compute the latitude of the
259 c - reference vector, in degrees, to send back as "rv0y".
260
261 c First, get the ratio of dx/dy in the lower left part of the
262 c - plot. Use 1% of the E/W span and apply that both ways.
263  ddum = 0.01d0*(be(i) - bw(i))
264  ydum = bs(i) + ddum
265  xdum = bw(i) + ddum
266
267  xw = bw(i) * d2r
268  xe = xdum * d2r
269
270  dx = (xe - xw)
271
272  q1 = dtan(ydum*d2r)
273  q2 = 1.d0 / dcos(ydum*d2r) ! Secant
274  yn = log(q1 + q2)
275
276  q1 = dtan(bs(i)*d2r)
277  q2 = 1.d0 / dcos(bs(i)*d2r) ! Secant
278  ys = log(q1 + q2)
279
280 c - dy is in radians
281  dy = yn - ys
282
283 c - Now I have a new "dx" and "dy" representing the lower left 1% of the plot,
284 c - in radians/radians. Compute the new ratio for this part.
285  ratio = dx/dy
286
287 c - Since radians to inches does NOT change in the E/W direction, and since
288 c - I know exactly how wide the plot is, I can compute a radians/inch conversion
289 c - in E/W, for 1% of the width of the plot.
290  r2iew = (0.01d0 * jm(i)) / dx
291
292 c - The ratio of dx/dy, is going to be identical to the ratio of r2iew/r2ins:
293  r2ins = r2iew / ratio
294
295 c - Convert r2ins into degrees-to-inches
296  d2ins = r2ins * d2r
297
298 c - Now, the south edge is known in degrees. I want the reference vector
299 c - to be 1/2 inch below that. I want the reference vector label to
300 c - be an additional 1/4 inches down.
301  rv0y(i) = bs(i) - (0.50d0)/d2ins
302  rl0y(i) = bs(i) - (0.65d0)/d2ins
303
304 c - And set the reference vector's left edge to equal the left edge of the plot
305  rv0x(i) = bw(i)
306
307 c - Flags for debugging
308  write(6,8) trim(fn(i))
309  8 format('***************************',/,
310  * 'Reference Vector Computations for',
311  * 'Sub-Region: ',a)
312
313  write(6,*) ' yn,ys,dy = ',yn,ys,dy
314  write(6,*) ' xe,xw,dx = ',xe,xw,dx
315  write(6,*) ' E/W width (jm) = ',jm(i)
316  write(6,*) ' SW ratio = ',ratio
317  write(6,*) ' SW bounds,Lat = ',ydum,bs(i)
318  write(6,*) ' SW bounds,Lon = ',xdum,bw(i)
319  write(6,*) ' r2i E/W = ',r2iew
320  write(6,*) ' r2i N/S = ',r2ins
321  write(6,*) ' d2i N/S = ',d2ins
322  write(6,*) ' Southwest corner = ',bs(i),bw(i)
323  write(6,*) ' Ref Vector = ',rv0y(i),rv0x(i)
324  write(6,*) ' Ref Label = ',rl0y(i),rv0x(i)
325
326  goto 1
327
328  2 nplots = i
329
330  close(ifile)
331
332  if(nplots.eq.0)then
333  write(6,100)trim(region)
334  stop 10001
335  endif
336  100 format(6x,'FATAL Subroutine getmapbounds: Unknown region: ',a)
337  return
338  end
subroutine getmapbounds(mapflag, maxplots, region, nplots, olddtm, newdtm, bw, be, bs, bn, jm, b1, b2, fn, lrv, rv0x, rv0y, rl0y)
Subroutine to collect up the MAP boundaries for use in creating NADCON 5.
Definition: getmapbounds.f:73