27 DOUBLE PRECISION :: latin1, latin2
28 DOUBLE PRECISION :: lad
29 DOUBLE PRECISION :: lon1, lat1
30 INTEGER :: projection_center_flag
31 END TYPE geo_proj_polar
34 DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
35 END TYPE geo_proj_rotated
37 TYPE geo_proj_stretched
38 DOUBLE PRECISION :: latitude_stretch_pole, longitude_stretch_pole, stretch_factor
39 END TYPE geo_proj_stretched
43 DOUBLE PRECISION :: rf, a
44 DOUBLE PRECISION :: f, e2, e1, ep2, e11, e12, e13, e14, e4, e6, ef0, ef1, ef2, ef3, k0
45 END TYPE geo_proj_ellips
48 CHARACTER(len=80) :: proj_type=cmiss
49 DOUBLE PRECISION :: xoff, yoff
50 DOUBLE PRECISION :: lov
51 TYPE(geo_proj_rotated
) :: rotated
52 TYPE(geo_proj_stretched
) :: stretched
53 TYPE(geo_proj_polar
) :: polar
55 TYPE(geo_proj_ellips
) :: ellips
61 MODULE PROCEDURE geo_proj_delete
66 MODULE PROCEDURE geo_proj_copy
71 MODULE PROCEDURE geo_proj_get_val
76 MODULE PROCEDURE geo_proj_set_val
81 MODULE PROCEDURE geo_proj_c_e
86 MODULE PROCEDURE geo_proj_write_unit
91 MODULE PROCEDURE geo_proj_read_unit
96 MODULE PROCEDURE geo_proj_display
102 MODULE PROCEDURE geo_proj_proj
108 MODULE PROCEDURE geo_proj_unproj
114 INTERFACE operator (==)
115 MODULE PROCEDURE geo_proj_eq
121 INTERFACE operator (/=)
122 MODULE PROCEDURE geo_proj_ne
126 INTEGER,
PARAMETER,
PUBLIC :: &
127 geo_proj_unit_degree = 0, & !< coordinate unit is degrees (longitude periodic over 360 degrees)
128 geo_proj_unit_meter = 1
131 INTEGER,
PARAMETER :: nellips = 41
136 INTEGER,
PARAMETER,
PUBLIC :: &
137 ellips_merit = 1, & !< constants for predefined ellipsoids MERIT 1983
138 ellips_sgs85 = 2, & !< Soviet Geodetic System 85
139 ellips_grs80 = 3, & !< GRS 1980(IUGG, 1980)
140 ellips_iau76 = 4, & !< IAU 1976
141 ellips_airy = 5, & !< Airy 1830
142 ellips_apl4_9 = 6, & !< Appl. Physics. 1965
143 ellips_nwl9d = 7, & !< Naval Weapons Lab., 1965
144 ellips_mod_airy = 8, & !< Modified Airy
145 ellips_andrae = 9, & !< Andrae 1876 (Den., Iclnd.)
146 ellips_aust_sa = 10, & !< Australian Natl & S. Amer. 1969
147 ellips_grs67 = 11, & !< GRS 67(IUGG 1967)
148 ellips_bessel = 12, & !< Bessel 1841
149 ellips_bess_nam = 13, & !< Bessel 1841 (Namibia)
150 ellips_clrk66 = 14, & !< Clarke 1866
151 ellips_clrk80 = 15, & !< Clarke 1880 mod.
152 ellips_cpm = 16, & !< Comm. des Poids et Mesures 1799
153 ellips_delmbr = 17, & !< Delambre 1810 (Belgium)
154 ellips_engelis = 18, & !< Engelis 1985
155 ellips_evrst30 = 19, & !< Everest 1830
156 ellips_evrst48 = 20, & !< Everest 1948
157 ellips_evrst56 = 21, & !< Everest 1956
158 ellips_evrst69 = 22, & !< Everest 1969
159 ellips_evrstss = 23, & !< Everest (Sabah & Sarawak)
160 ellips_fschr60 = 24, & !< Fischer (Mercury Datum) 1960
161 ellips_fschr60m = 25, & !< Modified Fischer 1960
162 ellips_fschr68 = 26, & !< Fischer 1968
163 ellips_helmert = 27, & !< Helmert 1906
164 ellips_hough = 28, & !< Hough
165 ellips_intl = 29, & !< International 1909 (Hayford)
166 ellips_krass = 30, & !< Krassovsky, 1942
167 ellips_kaula = 31, & !< Kaula 1961
168 ellips_lerch = 32, & !< Lerch 1979
169 ellips_mprts = 33, & !< Maupertius 1738
170 ellips_new_intl = 34, & !< New International 1967
171 ellips_plessis = 35, & !< Plessis 1817 (France)
172 ellips_seasia = 36, & !< Southeast Asia
173 ellips_walbeck = 37, & !< Walbeck
174 ellips_wgs60 = 38, & !< WGS 60
175 ellips_wgs66 = 39, & !< WGS 66
176 ellips_wgs72 = 40, & !< WGS 72
179 DOUBLE PRECISION,
PARAMETER,
PRIVATE :: &
222 DOUBLE PRECISION,
PARAMETER,
PRIVATE :: &
266 DOUBLE PRECISION,
PARAMETER,
PRIVATE :: k0=0.9996d0
269 PUBLIC geo_proj, geo_proj_rotated, geo_proj_stretched, geo_proj_polar, &
283 FUNCTION geo_proj_new(proj_type, lov, zone, xoff, yoff, &
284 longitude_south_pole, latitude_south_pole, angle_rotation, &
285 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
286 latin1, latin2, lad, projection_center_flag, &
287 ellips_smaj_axis, ellips_flatt, ellips_type) result(this)
288 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
289 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
290 INTEGER,
INTENT(in),
OPTIONAL :: zone
291 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
292 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
293 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
294 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
295 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
296 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
297 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
298 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
299 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
300 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
301 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
302 INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
303 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
304 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
305 INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
307 TYPE(geo_proj
) :: this
309 this%proj_type = optio_c(proj_type, len(this%proj_type))
312 this%lov = optio_d(lov)
313 CALL
set_val(this, zone=zone, lov=lov)
316 this%xoff = optio_d(xoff)
317 IF (.NOT.
c_e(this%xoff)) this%xoff = 0.0d0
318 this%yoff = optio_d(yoff)
319 IF (.NOT.
c_e(this%yoff)) this%yoff = 0.0d0
321 this%rotated%longitude_south_pole = optio_d(longitude_south_pole)
322 this%rotated%latitude_south_pole = optio_d(latitude_south_pole)
323 this%rotated%angle_rotation = optio_d(angle_rotation)
324 this%stretched%longitude_stretch_pole = optio_d(longitude_stretch_pole)
325 this%stretched%latitude_stretch_pole = optio_d(latitude_stretch_pole)
326 this%stretched%stretch_factor = optio_d(stretch_factor)
327 this%polar%latin1 = optio_d(latin1)
328 this%polar%latin2 = optio_d(latin2)
329 this%polar%lad = optio_d(lad)
330 this%polar%projection_center_flag = optio_l(projection_center_flag)
333 CALL ellips_compute(this%ellips)
334 CALL
set_val(this, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
335 ellips_type=ellips_type)
337 END FUNCTION geo_proj_new
342 SUBROUTINE ellips_compute(this, a, f)
343 TYPE(geo_proj_ellips
),
INTENT(inout) :: this
344 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: a, f
346 IF (present(a) .AND. present(f))
THEN
349 ELSE IF (present(a))
THEN
357 this%e2 = 2.0d0*this%f - this%f*this%f
358 this%e1 = this%f/(2.0d0 - this%f)
359 this%ep2 = this%e2/(1.0d0 - this%e2)
360 this%e11 = 3.0d0*this%e1/2.0d0 - 27.0d0*this%e1*this%e1*this%e1/32.0d0
361 this%e12 = 21.0d0*this%e1*this%e1/16.0d0 &
362 - 55.0d0*this%e1*this%e1*this%e1*this%e1/32.0d0
363 this%e13 = 151.0d0*this%e1*this%e1*this%e1/96.0d0
364 this%e14 = 1097.0d0*this%e1*this%e1*this%e1*this%e1/512.0d0
365 this%e4 = this%e2*this%e2
366 this%e6 = this%e2*this%e4
367 this%ef0 = this%a*(1.0d0 - 0.25d0*this%e2&
368 *(1.0d0+this%e2/16.0d0*(3.0d0 + 1.25d0*this%e2)))
369 this%ef1 = this%a*(0.375d0*this%e2 &
370 *(1.0d0 + 0.25d0*this%e2*(1.0d0 + 0.46875d0*this%e2)))
371 this%ef2 = this%a*(0.05859375d0*this%e2*this%e2*(1.0d0 + 0.75d0*this%e2))
372 this%ef3 = this%a*this%e2*this%e2*this%e2*35.0d0/3072.0d0
374 END SUBROUTINE ellips_compute
378 SUBROUTINE geo_proj_delete(this)
379 TYPE(geo_proj
),
INTENT(inout) :: this
381 this%proj_type = cmiss
385 this%rotated%longitude_south_pole = dmiss
386 this%rotated%latitude_south_pole = dmiss
387 this%rotated%angle_rotation = dmiss
388 this%stretched%longitude_stretch_pole = dmiss
389 this%stretched%latitude_stretch_pole = dmiss
390 this%stretched%stretch_factor = dmiss
391 this%polar%latin1 = dmiss
392 this%polar%latin2 = dmiss
393 this%polar%lad = dmiss
394 this%polar%projection_center_flag = imiss
396 END SUBROUTINE geo_proj_delete
400 SUBROUTINE geo_proj_copy(this, that)
401 TYPE(geo_proj
),
INTENT(in) :: this
402 TYPE(geo_proj
),
INTENT(out) :: that
406 END SUBROUTINE geo_proj_copy
409 SUBROUTINE geo_proj_get_val(this, &
410 proj_type, lov, zone, xoff, yoff, &
411 longitude_south_pole, latitude_south_pole, angle_rotation, &
412 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
413 latin1, latin2, lad, projection_center_flag, &
414 ellips_smaj_axis, ellips_flatt, ellips_type, unit)
415 TYPE(geo_proj
),
INTENT(in) :: this
416 CHARACTER(len=*),
OPTIONAL :: proj_type
417 DOUBLE PRECISION,
OPTIONAL :: lov
418 INTEGER,
OPTIONAL :: zone
419 DOUBLE PRECISION,
OPTIONAL :: xoff
420 DOUBLE PRECISION,
OPTIONAL :: yoff
421 DOUBLE PRECISION,
OPTIONAL :: longitude_south_pole
422 DOUBLE PRECISION,
OPTIONAL :: latitude_south_pole
423 DOUBLE PRECISION,
OPTIONAL :: angle_rotation
424 DOUBLE PRECISION,
OPTIONAL :: longitude_stretch_pole
425 DOUBLE PRECISION,
OPTIONAL :: latitude_stretch_pole
426 DOUBLE PRECISION,
OPTIONAL :: stretch_factor
427 DOUBLE PRECISION,
OPTIONAL :: latin1
428 DOUBLE PRECISION,
OPTIONAL :: latin2
429 DOUBLE PRECISION,
OPTIONAL :: lad
430 INTEGER,
OPTIONAL :: projection_center_flag
431 DOUBLE PRECISION,
OPTIONAL :: ellips_smaj_axis
432 DOUBLE PRECISION,
OPTIONAL :: ellips_flatt
433 INTEGER,
OPTIONAL :: ellips_type
434 INTEGER,
OPTIONAL :: unit
438 IF (present(proj_type)) proj_type = this%proj_type
440 IF (present(lov) .AND. present(zone))
THEN
441 zone = nint((this%lov + 183.0d0)/6.0d0)
442 lov = this%lov - zone*6.0d0 - 183.0d0
443 IF (
abs(lov) < 1.0d-6) lov = 0.0d-6
444 ELSE IF (present(lov))
THEN
446 ELSE IF (present(zone))
THEN
447 zone = nint((this%lov + 183.0d0)/6.0d0)
450 IF (present(xoff)) xoff = this%xoff
451 IF (present(yoff)) yoff = this%yoff
452 IF (present(longitude_south_pole)) longitude_south_pole = this%rotated%longitude_south_pole
453 IF (present(latitude_south_pole)) latitude_south_pole = this%rotated%latitude_south_pole
454 IF (present(angle_rotation)) angle_rotation = this%rotated%angle_rotation
455 IF (present(longitude_stretch_pole)) longitude_stretch_pole = this%stretched%longitude_stretch_pole
456 IF (present(latitude_stretch_pole)) latitude_stretch_pole = this%stretched%latitude_stretch_pole
457 IF (present(stretch_factor)) stretch_factor = this%stretched%stretch_factor
458 IF (present(latin1)) latin1 = this%polar%latin1
459 IF (present(latin2)) latin2 = this%polar%latin2
460 IF (present(lad)) lad = this%polar%lad
461 IF (present(projection_center_flag)) projection_center_flag = this%polar%projection_center_flag
463 IF (present(ellips_smaj_axis)) ellips_smaj_axis = this%ellips%a
464 IF (present(ellips_flatt)) ellips_flatt = this%ellips%f
465 IF (present(ellips_type))
THEN
468 IF (this%ellips%f == 1.0d0/rf(i) .AND. this%ellips%a == a(i))
THEN
475 IF (present(unit))
THEN
476 SELECT CASE(this%proj_type)
477 CASE(
"regular_ll",
"rotated_ll")
478 unit = geo_proj_unit_degree
479 CASE(
"lambert",
"polar_stereographic",
"UTM")
480 unit = geo_proj_unit_meter
487 END SUBROUTINE geo_proj_get_val
490 SUBROUTINE geo_proj_set_val(this, &
491 proj_type, lov, zone, xoff, yoff, &
492 longitude_south_pole, latitude_south_pole, angle_rotation, &
493 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
494 latin1, latin2, lad, projection_center_flag, &
495 ellips_smaj_axis, ellips_flatt, ellips_type)
496 TYPE(geo_proj
),
INTENT(inout) :: this
497 CHARACTER(len=*),
OPTIONAL :: proj_type
498 DOUBLE PRECISION,
OPTIONAL :: lov
499 INTEGER,
INTENT(in),
OPTIONAL :: zone
500 DOUBLE PRECISION,
OPTIONAL :: xoff
501 DOUBLE PRECISION,
OPTIONAL :: yoff
502 DOUBLE PRECISION,
OPTIONAL :: longitude_south_pole
503 DOUBLE PRECISION,
OPTIONAL :: latitude_south_pole
504 DOUBLE PRECISION,
OPTIONAL :: angle_rotation
505 DOUBLE PRECISION,
OPTIONAL :: longitude_stretch_pole
506 DOUBLE PRECISION,
OPTIONAL :: latitude_stretch_pole
507 DOUBLE PRECISION,
OPTIONAL :: stretch_factor
508 DOUBLE PRECISION,
OPTIONAL :: latin1
509 DOUBLE PRECISION,
OPTIONAL :: latin2
510 DOUBLE PRECISION,
OPTIONAL :: lad
511 INTEGER,
OPTIONAL :: projection_center_flag
512 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
513 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
514 INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
517 DOUBLE PRECISION :: llov
522 lzone = optio_i(zone)
523 IF (
c_e(llov) .AND.
c_e(lzone))
THEN
524 this%lov = llov + zone*6.0d0 - 183.0d0
525 ELSE IF (
c_e(llov))
THEN
527 ELSE IF (
c_e(lzone))
THEN
528 this%lov = lzone*6.0d0 - 183.0d0
532 IF (present(ellips_smaj_axis))
THEN
534 CALL ellips_compute(this%ellips, ellips_smaj_axis, ellips_flatt)
535 ELSE IF (present(ellips_type))
THEN
536 IF (ellips_type > 0 .AND. ellips_type <= nellips)
THEN
538 CALL ellips_compute(this%ellips, a(ellips_type), 1.0d0/rf(ellips_type))
540 CALL ellips_compute(this%ellips)
546 END SUBROUTINE geo_proj_set_val
549 FUNCTION geo_proj_c_e(this) RESULT(c_e)
550 TYPE(geo_proj
),
INTENT(in) :: this
553 c_e = this%proj_type /= cmiss
555 END FUNCTION geo_proj_c_e
562 SUBROUTINE geo_proj_read_unit(this, unit)
563 TYPE(geo_proj
),
INTENT(out) :: this
564 INTEGER,
INTENT(in) :: unit
566 CHARACTER(len=40) :: form
568 INQUIRE(unit, form=form)
569 IF (form ==
'FORMATTED')
THEN
575 END SUBROUTINE geo_proj_read_unit
582 SUBROUTINE geo_proj_write_unit(this, unit)
583 TYPE(geo_proj
),
INTENT(in) :: this
584 INTEGER,
INTENT(in) :: unit
586 CHARACTER(len=40) :: form
588 INQUIRE(unit, form=form)
589 IF (form ==
'FORMATTED')
THEN
595 END SUBROUTINE geo_proj_write_unit
598 SUBROUTINE geo_proj_display(this)
599 TYPE(geo_proj
),
INTENT(in) :: this
601 print*,
"<<<<<<<<<<<<<<< ",
t2c(this%proj_type,
"undefined projection"), &
604 IF (
c_e(this%xoff) .AND. this%xoff /= 0.0d0)
THEN
605 print*,
"False easting",this%xoff
607 IF (
c_e(this%yoff) .AND. this%yoff /= 0.0d0)
THEN
608 print*,
"False northing",this%yoff
611 IF (
c_e(this%lov))
THEN
612 print*,
"Central Meridian",this%lov
615 IF (this%proj_type ==
'rotated_ll' .OR. this%proj_type ==
'stretched_rotated_ll')
THEN
616 print*,
"Rotated projection:"
617 print*,
"lon of south pole",this%rotated%longitude_south_pole
618 print*,
"lat of south pole",this%rotated%latitude_south_pole
619 print*,
"angle of rotation",this%rotated%angle_rotation
622 IF (this%proj_type ==
'stretched_ll' .OR. this%proj_type ==
'stretched_rotated_ll')
THEN
623 print*,
"Stretched projection:"
624 print*,
"lon of stretch pole",this%stretched%longitude_stretch_pole
625 print*,
"lat of stretch pole",this%stretched%latitude_stretch_pole
626 print*,
"stretching factor",this%stretched%stretch_factor
629 IF (this%proj_type ==
'lambert' .OR. this%proj_type ==
'polar_stereographic')
THEN
630 print*,
"Polar projection:"
631 IF (
c_e(this%polar%latin1) .OR.
c_e(this%polar%latin2))
THEN
632 print*,
"lat of intersections",this%polar%latin1,this%polar%latin2
634 IF (
c_e(this%polar%lad))
THEN
635 print*,
"isometric latitude",this%polar%lad
637 IF (iand(this%polar%projection_center_flag, 128) == 0)
THEN
644 IF (this%ellips%f == 0.0d0)
THEN
645 print*,
"Spherical Earth:"
646 print*,
"Radius (m)",this%ellips%a
649 print*,
"Flattening",this%ellips%f
650 print*,
"Reverse of flattening",1.0d0/this%ellips%f
651 print*,
"Semi-major axis (m)",this%ellips%a
655 END SUBROUTINE geo_proj_display
660 ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
661 TYPE(geo_proj
),
INTENT(in) :: this
663 DOUBLE PRECISION,
INTENT(in) :: lon, lat
665 DOUBLE PRECISION,
INTENT(out) :: x, y
667 SELECT CASE(this%proj_type)
670 CALL proj_regular_ll(lon, lat, x, y)
673 CALL proj_rotated_ll(lon, lat, x, y, this%rotated%longitude_south_pole, &
674 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
677 CALL proj_lambert(lon, lat, x, y, this%polar%latin1, &
678 this%polar%latin2, this%lov, this%polar%lad, &
679 this%polar%projection_center_flag)
681 CASE(
"polar_stereographic")
682 CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
683 this%polar%lad, this%polar%projection_center_flag)
686 CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
694 END SUBROUTINE geo_proj_proj
699 ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
700 TYPE(geo_proj
),
INTENT(in) :: this
702 DOUBLE PRECISION,
INTENT(in) :: x, y
704 DOUBLE PRECISION,
INTENT(out) :: lon, lat
706 SELECT CASE(this%proj_type)
709 CALL unproj_regular_ll(x, y, lon, lat)
712 CALL unproj_rotated_ll(x, y, lon, lat, this%rotated%longitude_south_pole, &
713 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
716 CALL unproj_lambert(x, y, lon, lat, this%polar%latin1, &
717 this%polar%latin2, this%lov, this%polar%lad, &
718 this%polar%projection_center_flag)
720 CASE(
"polar_stereographic")
721 CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
722 this%polar%lad, this%polar%projection_center_flag)
725 CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
733 END SUBROUTINE geo_proj_unproj
736 ELEMENTAL FUNCTION geo_proj_eq(this, that) RESULT(eq)
737 TYPE(geo_proj
),
INTENT(in) :: this, that
740 eq = this%proj_type == that%proj_type .AND. this%xoff == that%xoff .AND. &
741 this%yoff == that%yoff .AND. geo_lon_eq(this%lov, that%lov) .AND. &
742 geo_lon_eq(this%rotated%longitude_south_pole, that%rotated%longitude_south_pole) .AND. &
743 this%rotated%latitude_south_pole == that%rotated%latitude_south_pole .AND. &
744 this%rotated%angle_rotation == that%rotated%angle_rotation .AND. &
745 this%stretched%latitude_stretch_pole == that%stretched%latitude_stretch_pole .AND. &
746 geo_lon_eq(this%stretched%longitude_stretch_pole, that%stretched%longitude_stretch_pole) .AND. &
747 this%stretched%stretch_factor == that%stretched%stretch_factor .AND. &
748 this%polar%latin1 == that%polar%latin1 .AND. &
749 this%polar%latin2 == that%polar%latin2 .AND. &
750 this%polar%lad == that%polar%lad .AND. &
751 this%polar%projection_center_flag == that%polar%projection_center_flag .AND. &
752 this%ellips%f == that%ellips%f .AND. this%ellips%a == that%ellips%a
754 END FUNCTION geo_proj_eq
757 ELEMENTAL FUNCTION geo_proj_ne(this, that) RESULT(ne)
758 TYPE(geo_proj
),
INTENT(in) :: this, that
761 ne = .NOT. (this == that)
763 END FUNCTION geo_proj_ne
768 ELEMENTAL FUNCTION geo_lon_eq(l1, l2) RESULT(eq)
769 DOUBLE PRECISION,
INTENT(in) :: l1
770 DOUBLE PRECISION,
INTENT(in) :: l2
775 eq = modulo(l2-l1, 360.0d0) == 0.0d0
777 END FUNCTION geo_lon_eq
783 ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
784 DOUBLE PRECISION,
INTENT(in) :: lon,lat
785 DOUBLE PRECISION,
INTENT(out) :: x,y
790 END SUBROUTINE proj_regular_ll
792 ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
793 DOUBLE PRECISION,
INTENT(in) :: x,y
794 DOUBLE PRECISION,
INTENT(out) :: lon,lat
799 END SUBROUTINE unproj_regular_ll
802 ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
803 longitude_south_pole, latitude_south_pole, angle_rotation)
804 DOUBLE PRECISION,
INTENT(in) :: lon,lat
805 DOUBLE PRECISION,
INTENT(out) :: x,y
806 DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
809 DOUBLE PRECISION :: cy0,sy0,rx,srx,crx,sy,cy,l_south_pole
813 l_south_pole = (latitude_south_pole+90.)*degrad
815 rx = degrad*(lon - longitude_south_pole)
819 sy0 = sin(l_south_pole)
820 cy0 = cos(l_south_pole)
825 x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
826 y = raddeg*asin(cy0*sy - sy0*cy*crx)
827 x = x + angle_rotation
829 END SUBROUTINE proj_rotated_ll
831 ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
832 longitude_south_pole, latitude_south_pole, angle_rotation)
833 DOUBLE PRECISION,
INTENT(in) :: x,y
834 DOUBLE PRECISION,
INTENT(out) :: lon,lat
835 DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
838 DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
840 xr = (x - angle_rotation)*degrad
843 l_south_pole = (latitude_south_pole+90.)*degrad
845 cy0 = cos(l_south_pole)
846 sy0 = sin(l_south_pole)
848 lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
849 lon = longitude_south_pole + &
850 raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
852 END SUBROUTINE unproj_rotated_ll
855 ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
856 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
857 DOUBLE PRECISION,
INTENT(in) :: lon,lat
858 DOUBLE PRECISION,
INTENT(out) :: x,y
859 DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
862 DOUBLE PRECISION :: csq
864 csq = stretch_factor**2
866 y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
867 (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
869 END SUBROUTINE proj_stretched_ll
871 ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
872 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
873 DOUBLE PRECISION,
INTENT(in) :: x,y
874 DOUBLE PRECISION,
INTENT(out) :: lon,lat
875 DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
878 DOUBLE PRECISION :: csq
880 csq = stretch_factor**2
883 lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
884 (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
886 END SUBROUTINE unproj_stretched_ll
897 ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
898 latin1, latin2, lov, lad, projection_center_flag)
899 DOUBLE PRECISION,
INTENT(in) :: lon,lat
900 DOUBLE PRECISION,
INTENT(out) :: x,y
901 DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
902 INTEGER,
INTENT(in) :: projection_center_flag
904 DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
905 DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
907 IF (iand(projection_center_flag, 128) == 0)
THEN
908 pollat = 90.d0*degrad
910 pollat = -90.d0*degrad
912 cs1 = cos(degrad*latin1)
913 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
915 IF (latin1 == latin2)
THEN
916 n = sin(degrad*latin1)
918 n = log(cs1/cos(degrad*latin2)) / &
919 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
921 f = cs1*cs2**n/n*rearth
922 angle = pi*.25d0 + pollat*.5d0
923 cot = cos(angle)/sin(angle)
930 angle = pi*.25d0 + degrad*lat*.5d0
931 cot = cos(angle)/sin(angle)
938 cs3 = degrad*n*(lon - lov)
941 y = ro0 - ro*cos(cs3)
943 END SUBROUTINE proj_lambert
945 ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
946 latin1, latin2, lov, lad, projection_center_flag)
947 DOUBLE PRECISION,
INTENT(in) :: x,y
948 DOUBLE PRECISION,
INTENT(out) :: lon,lat
949 DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
950 INTEGER,
INTENT(in) :: projection_center_flag
952 DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
953 DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
957 IF (iand(projection_center_flag, 128) == 0)
THEN
958 pollat = 90.d0*degrad
960 pollat = -90.d0*degrad
962 cs1 = cos(degrad*latin1)
963 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
965 IF (latin1 == latin2)
THEN
966 n = sin(degrad*latin1)
968 n = log(cs1/cos(degrad*latin2)) / &
969 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
971 f = cs1*cs2**n/n*rearth
972 angle = pi*.25d0 + pollat*.5d0
973 cot = cos(angle)/sin(angle)
980 ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n)
981 theta = raddeg*atan2(x, ro0-y)
984 lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
986 END SUBROUTINE unproj_lambert
990 ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
991 lov, lad, projection_center_flag)
992 DOUBLE PRECISION,
INTENT(in) :: lon,lat
993 DOUBLE PRECISION,
INTENT(out) :: x,y
994 DOUBLE PRECISION,
INTENT(in) :: lov, lad
995 INTEGER,
INTENT(in) :: projection_center_flag
997 DOUBLE PRECISION :: k, pollat
999 IF (iand(projection_center_flag, 128) == 0)
THEN
1000 pollat = 90.d0*degrad
1002 pollat = -90.d0*degrad
1005 k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1006 cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1007 x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1008 y = k*(cos(pollat)*sin(degrad*lat) - &
1009 sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1011 END SUBROUTINE proj_polar_stereographic
1013 ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1014 lov, lad, projection_center_flag)
1015 DOUBLE PRECISION,
INTENT(in) :: x,y
1016 DOUBLE PRECISION,
INTENT(out) :: lon,lat
1017 DOUBLE PRECISION,
INTENT(in) :: lov, lad
1018 INTEGER,
INTENT(in) :: projection_center_flag
1020 DOUBLE PRECISION :: ro, c, pollat
1022 IF (iand(projection_center_flag, 128) == 0)
THEN
1023 pollat = 90.d0*degrad
1025 pollat = -90.d0*degrad
1028 ro = sqrt(x**2 + y**2)
1029 c = 2.0d0*atan(ro/(2.0d0*rearth))
1030 lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1031 lon = lov + raddeg*atan2(x*sin(c), &
1032 (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1034 END SUBROUTINE unproj_polar_stereographic
1037 ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1038 DOUBLE PRECISION,
INTENT(in) :: lon,lat
1039 DOUBLE PRECISION,
INTENT(out) :: x,y
1040 DOUBLE PRECISION,
INTENT(in) :: lov
1041 DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1042 TYPE(geo_proj_ellips
),
INTENT(in) :: ellips
1044 DOUBLE PRECISION :: deltalon, p
1045 DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1048 deltalon = degrad*(lon - lov)
1056 n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1058 c = ellips%ep2*cosp*cosp
1063 m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1064 ellips%ef3*sin(6.0d0*p)
1074 x = k0*n*(a1 + (1.0d0 - t + c)*a3/6.0d0 &
1075 + (5.0d0 - 18.0d0*t + t2 + 72.0d0*c - 58.0d0*ellips%ep2) &
1076 *a5/120.0d0) + false_e
1077 y = k0*(m + n*tanp &
1078 *(a2/2.0d0 + (5.0d0 - t + 9.0d0*c + &
1079 4.0d0*c*c)*a4/24.0d0 + (61.0d0 - 58.0d0*t + t2 + &
1080 600.0d0*c - 330.0d0*ellips%ep2)*a6/720.0d0))
1083 END SUBROUTINE proj_utm
1085 ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1086 DOUBLE PRECISION,
INTENT(IN) :: x, y
1087 DOUBLE PRECISION,
INTENT(OUT) :: lon, lat
1088 DOUBLE PRECISION,
INTENT(in) :: lov
1089 DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1090 TYPE(geo_proj_ellips
),
INTENT(in) :: ellips
1092 DOUBLE PRECISION :: xm, ym, m, u, p1, c1, c2, t1, t2, n1, &
1093 sinp1, cosp1, tanp1, sin2p1, r0, r1, d, d2, d3, d4, d5, d6, p, l
1104 u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1105 3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1106 p1 = u + ellips%e11*sin(2.0d0*u) + ellips%e12*sin(4.0d0*u) + &
1107 ellips%e13*sin(6.0d0*u) + ellips%e14*sin(8.0d0*u)
1111 c1 = ellips%ep2*cosp1**2
1116 r0 = 1.0d0-ellips%e2*sin2p1
1117 n1 = ellips%a/sqrt(r0)
1119 r1 = r0/(1.0d0-ellips%e2)
1133 p = p1 - (r1*tanp1) * (d2/2.0d0 &
1134 - (5.0d0 + 3.0d0*t1 + 10.0d0*c1 - 4.0d0*c2 &
1135 - 9.0d0*ellips%ep2)*d4/24.0d0 &
1136 + (61.0d0 + 90.0d0*t1 + 298.0d0*c1 + 45.0d0*t2 &
1137 - 252d0*ellips%ep2 - 3.0d0*c2)*d6/720.0d0)
1139 l = (d - (1.0d0 + 2.0d0*t1 + c1)*d3/6.0d0 &
1140 + (5.0d0 - 2.0d0*c1 + 28.0d0*t1 - 3.0d0*c2 &
1141 + 8.0d0*ellips%ep2 + 24.0d0*t2)*d5/120.0d0)/cosp1
1142 lon = raddeg*l + lov
1144 END SUBROUTINE unproj_utm
1146 END MODULE geo_proj_class
Write the object on a formatted or unformatted file.
Read the object from a formatted or unformatted file.
Definitions of constants and functions for working with missing values.
Functions that return a trimmed CHARACTER representation of the input variable.
Compute forward coordinate transformation from geographical system to projected system.
Operatore di valore assoluto di un intervallo.
Compute backward coordinate transformation from projected system to geographical system.
Copy an object, creating a fully new instance.
Print a brief description on stdout.
Method for returning the contents of the object.
Destructors of the corresponding objects.
Utilities for CHARACTER variables.
Method for testing the existence of the object.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Method for setting the contents of the object.
Costanti fisiche (DOUBLEPRECISION).