-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathproject.F90
308 lines (275 loc) · 11.8 KB
/
project.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
program project
! This module requires Fortran 2003 or later
use iso_fortran_env, only : &
input_unit,output_unit,error_unit
use projection, only : &
PJ_ilatlonflag,PJ_iprojflag,PJ_k0,PJ_lam0,PJ_phi0,PJ_phi1,PJ_phi2,PJ_Re,&
PJ_Set_Proj_Params, &
PJ_proj_for, &
PJ_proj_inv
implicit none
integer :: nargs
integer :: iostatus
integer :: inlen
character(len=120) :: iomessage
character(len=100) :: arg
real(kind=8) :: inx ! input lon or x to convert
real(kind=8) :: iny ! input lat or y to convert
real(kind=8) :: outx ! output lon or x returned
real(kind=8) :: outy ! output lat or y returned
integer :: ProjDir
! CONUS 3.0-km Lambert Conformal
! 0 4 262.5 38.5 38.5 38.5 6371.229 #Proj flags and params
! proj +proj=lcc +lon_0=262.5 +lat_0=38.5 +lat_1=38.5 +lat_2=38.5 +R=6371.229
! CONUS 40-km Lambert Conformal and
! CONUS 12-km Lambert Conformal and
! CONUS 5.079-km Lambert Conformal
! 0 4 265.0 25.0 25.0 25.0 6371.229 #Proj flags and params
! proj +proj=lcc +lon_0=265.0 +lat_0=25.0 +lat_1=25.0 +lat_2=25.0 +R=6371.229
! NAM 32-km Lambert Conformal used by NARR (used PJ_Re=6367.470, not 6371.229)
! 0 4 -107.0 50.0 50.0 50.0 6367.47 #Proj flags and params
! proj +proj=lcc +lon_0=-107.0 +lat_0=50.0 +lat_1=50.0 +lat_2=50.0 +R=6367.47
! NAM 32-km Lambert Conformal
! HI 2.5-km Mercator
! 0 5 198.475 20.0 0.933 6371.229 #Proj flags and params
! proj +proj=merc +lat_ts=20.0 +lon_0=198.475 +R=6371.229
! NAM 6-km Polar Stereographic
! 0 1 -150.0 90.0 0.933 6371.229 #Proj flags and params
! proj +proj=stere +lon_0=210 +lat_0=90 +k_0=0.933 +R=6371.229
! NAM 3-km Polar Stereographic and
! NAM 90-km Polar Stereographic
! 0 1 -150.0 90.0 0.933 6371.229 #Proj flags and params
! proj +proj=stere +lon_0=210 +lat_0=90 +k_0=0.933 +R=6371.229
! NAM 45-km Polar Stereographic and
! NAM 11.25-km Polar Stereographic
! 0 1 -135.0 90.0 0.933 6371.229 #Proj flags and params
! proj +proj=stere +lon_0=225 +lat_0=90 +k_0=0.933 +R=6371.229
ProjDir = 0
#ifdef FORWARD
ProjDir = 1
#endif
#ifdef INVERSE
ProjDir = -1
#endif
if (ProjDir.eq.0) then
write(error_unit,*)"project.F90 must be compiled with either -DFORWARD or -DINVERSE"
write(error_unit,*)"Forward or inverse not set"
stop 1
endif
! TEST READ COMMAND LINE ARGUMENTS
nargs = command_argument_count()
if (nargs.lt.4) then
write(error_unit,*)"ERROR: Too few arguments. Enter lon lat IsLatLon ProjFlag ..."
stop 1
endif
call get_command_argument(number=1, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)inx
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (1)"
write(error_unit,*)" inx = ",inx
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=2, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)iny
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (2)"
write(error_unit,*)" iny = ",iny
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=3, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_ilatlonflag
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (3)"
write(error_unit,*)" PJ_ilatlonflag = ",PJ_ilatlonflag
write(error_unit,*)iomessage
stop 1
endif
if (PJ_ilatlonflag.eq.1) then
! coordinates are in lon/lat
stop 0
elseif (PJ_ilatlonflag.eq.0) then
! coordinates are projected, read the projection flag
call get_command_argument(number=4, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_iprojflag
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (4)"
write(error_unit,*)" PJ_iprojflag = ",PJ_iprojflag
write(error_unit,*)iomessage
stop 1
endif
if (PJ_iprojflag.ne.0.and.PJ_iprojflag.ne.1.and. &
PJ_iprojflag.ne.2.and.PJ_iprojflag.ne.3.and. &
PJ_iprojflag.ne.4.and.PJ_iprojflag.ne.5) then
write(error_unit,*)"Unrecognized projection flag"
stop 1
endif
else
! PJ_ilatlonflag is not 0 or 1, stopping program
write(error_unit,*)"Unrecognized latlonflag"
stop 1
endif
select case (PJ_iprojflag)
case(0)
! Non-geographic projection, (x,y) only
PJ_k0 = 0.0_8
PJ_Re = 6371.229_8
PJ_lam0 = 0.0_8
PJ_phi0 = 0.0_8
PJ_phi1 = 0.0_8
PJ_phi2 = 0.0_8
write(output_unit,*)"Both PJ_ilatlonflag and PJ_iprojflag are 0"
write(output_unit,*)"No geographic projection used"
case(1)
! Polar stereographic
!read(linebuffer,*)PJ_ilatlonflag,PJ_iprojflag,PJ_lam0,PJ_phi0,PJ_k0,PJ_Re
if (nargs.lt.8) then
write(error_unit,*)"Enter lon lat IsLatLon ProjFlag lam0 phi0 k0 radius"
stop 1
endif
call get_command_argument(number=5, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_lam0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (5)"
write(error_unit,*)" PJ_lam0 = ",PJ_lam0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=6, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_phi0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (6)"
write(error_unit,*)" PJ_phi0 = ",PJ_phi0
write(error_unit,*)iomessage
stop 1
endif
PJ_phi1 = PJ_phi0 ! Set the truescale lat to be phi0 with scale
! determined by k0
call get_command_argument(number=7, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_k0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (7)"
write(error_unit,*)" PJ_k0 = ",PJ_k0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=8, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_Re
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (8)"
write(error_unit,*)" PJ_Re = ",PJ_Re
write(error_unit,*)iomessage
stop 1
endif
case(2)
! Albers Equal Area
write(output_unit,*)"WARNING: Albers not yet verified"
!read(linebuffer,*)PJ_ilatlonflag,PJ_iprojflag,PJ_lam0,PJ_phi0,PJ_phi1,PJ_phi2
if (nargs.lt.8) then
write(error_unit,*)"Enter lon lat IsLatLon ProjFlag lam0 phi0 phi1 phi2"
stop 1
endif
case(3)
! UTM
write(output_unit,*)"WARNING: UTM not yet verified"
!read(linebuffer,*)PJ_ilatlonflag,PJ_iprojflag,izone,inorth
if (nargs.lt.6) then
write(error_unit,*)"Enter lon lat IsLatLon ProjFlag izone inorth"
stop 1
endif
case(4)
! Lambert conformal conic (NARR, NAM218, NAM221)
!read(linebuffer,*)PJ_ilatlonflag,PJ_iprojflag,PJ_lam0, &
! PJ_phi0,PJ_phi1,PJ_phi2,PJ_Re
if (nargs.lt.9) then
write(error_unit,*)"Enter lon lat IsLatLon ProjFlag lam0 phi0 phi1 phi2 radius"
stop 1
endif
call get_command_argument(number=5, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_lam0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (5)"
write(error_unit,*)" PJ_lam0 = ",PJ_lam0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=6, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_phi0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (6)"
write(error_unit,*)" PJ_phi0 = ",PJ_phi0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=7, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_phi1
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (7)"
write(error_unit,*)" PJ_phi1 = ",PJ_phi1
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=8, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_phi2
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (8)"
write(error_unit,*)" PJ_phi2 = ",PJ_phi2
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=9, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_Re
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (9)"
write(error_unit,*)" PJ_Re = ",PJ_Re
write(error_unit,*)iomessage
stop 1
endif
case(5)
! Mercator (NAM196)
!read(linebuffer,*)PJ_ilatlonflag,PJ_iprojflag,PJ_lam0,PJ_phi0,PJ_Rd
if (nargs.lt.7) then
write(error_unit,*)"Enter lon lat IsLatLon ProjFlag lam0 phi0 radius"
stop 1
endif
call get_command_argument(number=5, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_lam0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (5)"
write(error_unit,*)" PJ_lam0 = ",PJ_lam0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=6, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_phi0
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (6)"
write(error_unit,*)" PJ_phi0 = ",PJ_phi0
write(error_unit,*)iomessage
stop 1
endif
call get_command_argument(number=7, value=arg, length=inlen, status=iostatus)
read(arg,*,iostat=iostatus,iomsg=iomessage)PJ_Re
if (iostatus.ne.0) then
write(error_unit,*)"ERROR: could not read command-line argument (7)"
write(error_unit,*)" PJ_Re = ",PJ_Re
write(error_unit,*)iomessage
stop 1
endif
end select
if (ProjDir.eq.1) then
! Forward Projection: assuming inx,iny are lon,lat
call PJ_proj_for(inx,iny, &
PJ_iprojflag,&
PJ_lam0,PJ_phi0,PJ_phi1,PJ_phi2,PJ_k0,PJ_Re,&
outx,outy)
else
! Inverse Projection: assuming inx,iny are x,y
call PJ_proj_inv(inx,iny, &
PJ_iprojflag,&
PJ_lam0,PJ_phi0,PJ_phi1,PJ_phi2,PJ_k0,PJ_Re,&
outx,outy)
endif
write(output_unit,*)real(outx,kind=4),real(outy,kind=4)
end program project