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%proj_type ==
'mercator')
THEN
645 IF (
c_e(this%polar%lad))
THEN
646 print*,
"isometric latitude",this%polar%lad
650 IF (this%ellips%f == 0.0d0)
THEN
651 print*,
"Spherical Earth:"
652 print*,
"Radius (m)",this%ellips%a
655 print*,
"Flattening",this%ellips%f
656 print*,
"Reverse of flattening",1.0d0/this%ellips%f
657 print*,
"Semi-major axis (m)",this%ellips%a
661 END SUBROUTINE geo_proj_display
666 ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
667 TYPE(geo_proj
),
INTENT(in) :: this
669 DOUBLE PRECISION,
INTENT(in) :: lon, lat
671 DOUBLE PRECISION,
INTENT(out) :: x, y
673 SELECT CASE(this%proj_type)
676 CALL proj_regular_ll(lon, lat, x, y)
679 CALL proj_rotated_ll(lon, lat, x, y, this%rotated%longitude_south_pole, &
680 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
683 CALL proj_lambert(lon, lat, x, y, this%polar%latin1, &
684 this%polar%latin2, this%lov, this%polar%lad, &
685 this%polar%projection_center_flag)
687 CASE(
"polar_stereographic")
688 CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
689 this%polar%lad, this%polar%projection_center_flag)
692 CALL proj_mercator(lon, lat, x, y, this%lov, this%polar%lad)
695 CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
703 END SUBROUTINE geo_proj_proj
708 ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
709 TYPE(geo_proj
),
INTENT(in) :: this
711 DOUBLE PRECISION,
INTENT(in) :: x, y
713 DOUBLE PRECISION,
INTENT(out) :: lon, lat
715 SELECT CASE(this%proj_type)
718 CALL unproj_regular_ll(x, y, lon, lat)
721 CALL unproj_rotated_ll(x, y, lon, lat, this%rotated%longitude_south_pole, &
722 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
725 CALL unproj_lambert(x, y, lon, lat, this%polar%latin1, &
726 this%polar%latin2, this%lov, this%polar%lad, &
727 this%polar%projection_center_flag)
729 CASE(
"polar_stereographic")
730 CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
731 this%polar%lad, this%polar%projection_center_flag)
734 CALL unproj_mercator(x, y, lon, lat, this%lov, this%polar%lad)
737 CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
745 END SUBROUTINE geo_proj_unproj
748 ELEMENTAL FUNCTION geo_proj_eq(this, that) RESULT(eq)
749 TYPE(geo_proj
),
INTENT(in) :: this, that
752 eq = this%proj_type == that%proj_type .AND. this%xoff == that%xoff .AND. &
753 this%yoff == that%yoff .AND. geo_lon_eq(this%lov, that%lov) .AND. &
754 geo_lon_eq(this%rotated%longitude_south_pole, that%rotated%longitude_south_pole) .AND. &
755 this%rotated%latitude_south_pole == that%rotated%latitude_south_pole .AND. &
756 this%rotated%angle_rotation == that%rotated%angle_rotation .AND. &
757 this%stretched%latitude_stretch_pole == that%stretched%latitude_stretch_pole .AND. &
758 geo_lon_eq(this%stretched%longitude_stretch_pole, that%stretched%longitude_stretch_pole) .AND. &
759 this%stretched%stretch_factor == that%stretched%stretch_factor .AND. &
760 this%polar%latin1 == that%polar%latin1 .AND. &
761 this%polar%latin2 == that%polar%latin2 .AND. &
762 this%polar%lad == that%polar%lad .AND. &
763 this%polar%projection_center_flag == that%polar%projection_center_flag .AND. &
764 this%ellips%f == that%ellips%f .AND. this%ellips%a == that%ellips%a
766 END FUNCTION geo_proj_eq
769 ELEMENTAL FUNCTION geo_proj_ne(this, that) RESULT(ne)
770 TYPE(geo_proj
),
INTENT(in) :: this, that
773 ne = .NOT. (this == that)
775 END FUNCTION geo_proj_ne
780 ELEMENTAL FUNCTION geo_lon_eq(l1, l2) RESULT(eq)
781 DOUBLE PRECISION,
INTENT(in) :: l1
782 DOUBLE PRECISION,
INTENT(in) :: l2
787 eq = modulo(l2-l1, 360.0d0) == 0.0d0
789 END FUNCTION geo_lon_eq
795 ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
796 DOUBLE PRECISION,
INTENT(in) :: lon,lat
797 DOUBLE PRECISION,
INTENT(out) :: x,y
802 END SUBROUTINE proj_regular_ll
804 ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
805 DOUBLE PRECISION,
INTENT(in) :: x,y
806 DOUBLE PRECISION,
INTENT(out) :: lon,lat
811 END SUBROUTINE unproj_regular_ll
814 ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
815 longitude_south_pole, latitude_south_pole, angle_rotation)
816 DOUBLE PRECISION,
INTENT(in) :: lon,lat
817 DOUBLE PRECISION,
INTENT(out) :: x,y
818 DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
821 DOUBLE PRECISION :: cy0,sy0,rx,srx,crx,sy,cy,l_south_pole
825 l_south_pole = (latitude_south_pole+90.)*degrad
827 rx = degrad*(lon - longitude_south_pole)
831 sy0 = sin(l_south_pole)
832 cy0 = cos(l_south_pole)
837 x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
838 y = raddeg*asin(cy0*sy - sy0*cy*crx)
839 x = x + angle_rotation
841 END SUBROUTINE proj_rotated_ll
843 ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
844 longitude_south_pole, latitude_south_pole, angle_rotation)
845 DOUBLE PRECISION,
INTENT(in) :: x,y
846 DOUBLE PRECISION,
INTENT(out) :: lon,lat
847 DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
850 DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
852 xr = (x - angle_rotation)*degrad
855 l_south_pole = (latitude_south_pole+90.)*degrad
857 cy0 = cos(l_south_pole)
858 sy0 = sin(l_south_pole)
860 lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
861 lon = longitude_south_pole + &
862 raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
864 END SUBROUTINE unproj_rotated_ll
867 ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
868 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
869 DOUBLE PRECISION,
INTENT(in) :: lon,lat
870 DOUBLE PRECISION,
INTENT(out) :: x,y
871 DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
874 DOUBLE PRECISION :: csq
876 csq = stretch_factor**2
878 y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
879 (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
881 END SUBROUTINE proj_stretched_ll
883 ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
884 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
885 DOUBLE PRECISION,
INTENT(in) :: x,y
886 DOUBLE PRECISION,
INTENT(out) :: lon,lat
887 DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
890 DOUBLE PRECISION :: csq
892 csq = stretch_factor**2
895 lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
896 (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
898 END SUBROUTINE unproj_stretched_ll
909 ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
910 latin1, latin2, lov, lad, projection_center_flag)
911 DOUBLE PRECISION,
INTENT(in) :: lon,lat
912 DOUBLE PRECISION,
INTENT(out) :: x,y
913 DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
914 INTEGER,
INTENT(in) :: projection_center_flag
916 DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
917 DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
919 IF (iand(projection_center_flag, 128) == 0)
THEN
920 pollat = 90.d0*degrad
922 pollat = -90.d0*degrad
924 cs1 = cos(degrad*latin1)
925 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
927 IF (latin1 == latin2)
THEN
928 n = sin(degrad*latin1)
930 n = log(cs1/cos(degrad*latin2)) / &
931 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
933 f = cs1*cs2**n/n*rearth
934 angle = pi*.25d0 + pollat*.5d0
935 cot = cos(angle)/sin(angle)
942 angle = pi*.25d0 + degrad*lat*.5d0
943 cot = cos(angle)/sin(angle)
950 cs3 = degrad*n*(lon - lov)
953 y = ro0 - ro*cos(cs3)
955 END SUBROUTINE proj_lambert
957 ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
958 latin1, latin2, lov, lad, projection_center_flag)
959 DOUBLE PRECISION,
INTENT(in) :: x,y
960 DOUBLE PRECISION,
INTENT(out) :: lon,lat
961 DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
962 INTEGER,
INTENT(in) :: projection_center_flag
964 DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
965 DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
969 IF (iand(projection_center_flag, 128) == 0)
THEN
970 pollat = 90.d0*degrad
972 pollat = -90.d0*degrad
974 cs1 = cos(degrad*latin1)
975 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
977 IF (latin1 == latin2)
THEN
978 n = sin(degrad*latin1)
980 n = log(cs1/cos(degrad*latin2)) / &
981 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
983 f = cs1*cs2**n/n*rearth
984 angle = pi*.25d0 + pollat*.5d0
985 cot = cos(angle)/sin(angle)
992 ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n)
993 theta = raddeg*atan2(x, ro0-y)
996 lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
998 END SUBROUTINE unproj_lambert
1002 ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
1003 lov, lad, projection_center_flag)
1004 DOUBLE PRECISION,
INTENT(in) :: lon,lat
1005 DOUBLE PRECISION,
INTENT(out) :: x,y
1006 DOUBLE PRECISION,
INTENT(in) :: lov, lad
1007 INTEGER,
INTENT(in) :: projection_center_flag
1009 DOUBLE PRECISION :: k, pollat
1011 IF (iand(projection_center_flag, 128) == 0)
THEN
1012 pollat = 90.d0*degrad
1014 pollat = -90.d0*degrad
1017 k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1018 cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1019 x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1020 y = k*(cos(pollat)*sin(degrad*lat) - &
1021 sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1023 END SUBROUTINE proj_polar_stereographic
1025 ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1026 lov, lad, projection_center_flag)
1027 DOUBLE PRECISION,
INTENT(in) :: x,y
1028 DOUBLE PRECISION,
INTENT(out) :: lon,lat
1029 DOUBLE PRECISION,
INTENT(in) :: lov, lad
1030 INTEGER,
INTENT(in) :: projection_center_flag
1032 DOUBLE PRECISION :: ro, c, pollat
1034 IF (iand(projection_center_flag, 128) == 0)
THEN
1035 pollat = 90.d0*degrad
1037 pollat = -90.d0*degrad
1040 ro = sqrt(x**2 + y**2)
1041 c = 2.0d0*atan(ro/(2.0d0*rearth))
1042 lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1043 lon = lov + raddeg*atan2(x*sin(c), &
1044 (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1046 END SUBROUTINE unproj_polar_stereographic
1050 ELEMENTAL SUBROUTINE proj_mercator(lon, lat, x, y, lov, lad)
1051 DOUBLE PRECISION,
INTENT(in) :: lon,lat
1052 DOUBLE PRECISION,
INTENT(out) :: x,y
1053 DOUBLE PRECISION,
INTENT(in) :: lov
1054 DOUBLE PRECISION,
INTENT(in) :: lad
1056 DOUBLE PRECISION :: scx, scy
1058 scy = cos(degrad*lad)*rearth
1061 y = log(tan(0.25d0*pi + 0.5d0*lat*degrad))*scy
1063 END SUBROUTINE proj_mercator
1065 ELEMENTAL SUBROUTINE unproj_mercator(x, y, lon, lat, lov, lad)
1066 DOUBLE PRECISION,
INTENT(in) :: x,y
1067 DOUBLE PRECISION,
INTENT(out) :: lon,lat
1068 DOUBLE PRECISION,
INTENT(in) :: lov
1069 DOUBLE PRECISION,
INTENT(in) :: lad
1071 DOUBLE PRECISION :: scx, scy
1073 scy = cos(degrad*lad)*rearth
1076 lat = 2.0d0*atan(exp(y/scy))*raddeg-90.0d0
1078 END SUBROUTINE unproj_mercator
1081 ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1082 DOUBLE PRECISION,
INTENT(in) :: lon,lat
1083 DOUBLE PRECISION,
INTENT(out) :: x,y
1084 DOUBLE PRECISION,
INTENT(in) :: lov
1085 DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1086 TYPE(geo_proj_ellips
),
INTENT(in) :: ellips
1088 DOUBLE PRECISION :: deltalon, p
1089 DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1092 deltalon = degrad*(lon - lov)
1100 n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1102 c = ellips%ep2*cosp*cosp
1107 m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1108 ellips%ef3*sin(6.0d0*p)
1118 x = k0*n*(a1 + (1.0d0 - t + c)*a3/6.0d0 &
1119 + (5.0d0 - 18.0d0*t + t2 + 72.0d0*c - 58.0d0*ellips%ep2) &
1120 *a5/120.0d0) + false_e
1121 y = k0*(m + n*tanp &
1122 *(a2/2.0d0 + (5.0d0 - t + 9.0d0*c + &
1123 4.0d0*c*c)*a4/24.0d0 + (61.0d0 - 58.0d0*t + t2 + &
1124 600.0d0*c - 330.0d0*ellips%ep2)*a6/720.0d0))
1127 END SUBROUTINE proj_utm
1129 ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1130 DOUBLE PRECISION,
INTENT(IN) :: x, y
1131 DOUBLE PRECISION,
INTENT(OUT) :: lon, lat
1132 DOUBLE PRECISION,
INTENT(in) :: lov
1133 DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1134 TYPE(geo_proj_ellips
),
INTENT(in) :: ellips
1136 DOUBLE PRECISION :: xm, ym, m, u, p1, c1, c2, t1, t2, n1, &
1137 sinp1, cosp1, tanp1, sin2p1, r0, r1, d, d2, d3, d4, d5, d6, p, l
1148 u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1149 3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1150 p1 = u + ellips%e11*sin(2.0d0*u) + ellips%e12*sin(4.0d0*u) + &
1151 ellips%e13*sin(6.0d0*u) + ellips%e14*sin(8.0d0*u)
1155 c1 = ellips%ep2*cosp1**2
1160 r0 = 1.0d0-ellips%e2*sin2p1
1161 n1 = ellips%a/sqrt(r0)
1163 r1 = r0/(1.0d0-ellips%e2)
1177 p = p1 - (r1*tanp1) * (d2/2.0d0 &
1178 - (5.0d0 + 3.0d0*t1 + 10.0d0*c1 - 4.0d0*c2 &
1179 - 9.0d0*ellips%ep2)*d4/24.0d0 &
1180 + (61.0d0 + 90.0d0*t1 + 298.0d0*c1 + 45.0d0*t2 &
1181 - 252d0*ellips%ep2 - 3.0d0*c2)*d6/720.0d0)
1183 l = (d - (1.0d0 + 2.0d0*t1 + c1)*d3/6.0d0 &
1184 + (5.0d0 - 2.0d0*c1 + 28.0d0*t1 - 3.0d0*c2 &
1185 + 8.0d0*ellips%ep2 + 24.0d0*t2)*d5/120.0d0)/cosp1
1186 lon = raddeg*l + lov
1188 END SUBROUTINE unproj_utm
1190 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).