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
Write the object on a formatted or unformatted file.
Read the object from a formatted or unformatted file.
Compute forward coordinate transformation from geographical system to projected system.
Set of functions that return a trimmed CHARACTER representation of the input variable.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
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.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Destructors of the corresponding objects.
Method for testing the existence of the object.
Utilities for CHARACTER variables.
Method for setting the contents of the object.