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
144 ellips_mod_airy = 8, &
146 ellips_aust_sa = 10, &
148 ellips_bessel = 12, &
149 ellips_bess_nam = 13, &
150 ellips_clrk66 = 14, &
151 ellips_clrk80 = 15, &
153 ellips_delmbr = 17, &
154 ellips_engelis = 18, &
155 ellips_evrst30 = 19, &
156 ellips_evrst48 = 20, &
157 ellips_evrst56 = 21, &
158 ellips_evrst69 = 22, &
159 ellips_evrstss = 23, &
160 ellips_fschr60 = 24, &
161 ellips_fschr60m = 25, &
162 ellips_fschr68 = 26, &
163 ellips_helmert = 27, &
170 ellips_new_intl = 34, &
171 ellips_plessis = 35, &
172 ellips_seasia = 36, &
173 ellips_walbeck = 37, &
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
Set of functions that return a trimmed CHARACTER representation of the input variable.
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Destructors of the corresponding objects.
Print a brief description on stdout.
Method for returning the contents of the object.
Compute forward coordinate transformation from geographical system to projected system.
Read the object from a formatted or unformatted file.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Write the object on a formatted or unformatted file.
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.