com/gradoservice/proj4as/proj/ProjLaea.as

   1 /*******************************************************************************\r
   2 NAME                  LAMBERT AZIMUTHAL EQUAL-AREA\r
   3  \r
   4 PURPOSE:        Transforms input longitude and latitude to Easting and\r
   5                 Northing for the Lambert Azimuthal Equal-Area projection.  The\r
   6                 longitude and latitude must be in radians.  The Easting\r
   7                 and Northing values will be returned in meters.\r
   8 \r
   9 PROGRAMMER              DATE            \r
  10 ----------              ----           \r
  11 D. Steinwand, EROS      March, 1991   \r
  12 \r
  13 This function was adapted from the Lambert Azimuthal Equal Area projection\r
  14 code (FORTRAN) in the General Cartographic Transformation Package software\r
  15 which is available from the U.S. Geological Survey National Mapping Division.\r
  16  \r
  17 ALGORITHM REFERENCES\r
  18 \r
  19 1.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,\r
  20     The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.\r
  21 \r
  22 2.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological\r
  23     Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United\r
  24     State Government Printing Office, Washington D.C., 1987.\r
  25 \r
  26 3.  "Software Documentation for GCTP General Cartographic Transformation\r
  27     Package", U.S. Geological Survey National Mapping Division, May 1982.\r
  28 *******************************************************************************/\r
  29 \r
  30 package com.gradoservice.proj4as.proj\r
  31 {\r
  32         import com.gradoservice.proj4as.ProjPoint;\r
  33         import com.gradoservice.proj4as.ProjConstants;\r
  34         import com.gradoservice.proj4as.Datum;\r
  35                 \r
  36         public class ProjLaea extends AbstractProjProjection\r
  37         {\r
  38                 private var sin_lat_o:Number;\r
  39                 private var cos_lat_o:Number;\r
  40                 \r
  41                 public function ProjLaea(data:ProjParams)\r
  42                 {\r
  43                         super(data);\r
  44                 }\r
  45                 \r
  46                 /* Initialize the Lambert Azimuthal Equal Area projection\r
  47                   ------------------------------------------------------*/\r
  48                   override public function init():void\r
  49                   {\r
  50                     this.sin_lat_o=Math.sin(this.lat0);\r
  51                     this.cos_lat_o=Math.cos(this.lat0);\r
  52                   }\r
  53                 \r
  54                 /* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y\r
  55                   -----------------------------------------------------------------------*/\r
  56                   override public function forward(p:ProjPoint):ProjPoint\r
  57                   {\r
  58                 \r
  59                     /* Forward equations\r
  60                       -----------------*/\r
  61                     var lon:Number=p.x;\r
  62                     var lat:Number=p.y;\r
  63                     var delta_lon:Number = ProjConstants.adjust_lon(lon - this.long0);\r
  64                 \r
  65                     //v 1.0\r
  66                     var sin_lat:Number=Math.sin(lat);\r
  67                     var cos_lat:Number=Math.cos(lat);\r
  68                 \r
  69                     var sin_delta_lon:Number=Math.sin(delta_lon);\r
  70                     var cos_delta_lon:Number=Math.cos(delta_lon);\r
  71                 \r
  72                     var g:Number =this.sin_lat_o * sin_lat +this.cos_lat_o * cos_lat * cos_delta_lon;\r
  73                     if (g == -1.0) {\r
  74                       trace("laea:fwd:Point projects to a circle of radius ");\r
  75                       return null;\r
  76                     }\r
  77                     var ksp:Number = this.a * Math.sqrt(2.0 / (1.0 + g));\r
  78                     var x:Number = ksp * cos_lat * sin_delta_lon + this.x0;\r
  79                     var y:Number = ksp * (this.cos_lat_o * sin_lat - this.sin_lat_o * cos_lat * cos_delta_lon) + this.y0;\r
  80                     p.x = x;\r
  81                     p.y = y\r
  82                     return p;\r
  83                   }//lamazFwd()\r
  84                 \r
  85                 /* Inverse equations\r
  86                   -----------------*/\r
  87                   override public function inverse(p:ProjPoint):ProjPoint\r
  88                   {\r
  89                     p.x -= this.x0;\r
  90                     p.y -= this.y0;\r
  91                 \r
  92                     var Rh:Number = Math.sqrt(p.x *p.x +p.y * p.y);\r
  93                     var temp:Number = Rh / (2.0 * this.a);\r
  94                 \r
  95                     if (temp > 1) {\r
  96                       trace("laea:Inv:DataError");\r
  97                       return null;\r
  98                     }\r
  99                 \r
 100                     var z:Number = 2.0 * ProjConstants.asinz(temp);\r
 101                     var sin_z:Number=Math.sin(z);\r
 102                     var cos_z:Number=Math.cos(z);\r
 103                 \r
 104                     var lon:Number =this.long0;\r
 105                     if (Math.abs(Rh) > ProjConstants.EPSLN) {\r
 106                        var lat:Number = ProjConstants.asinz(this.sin_lat_o * cos_z +this. cos_lat_o * sin_z *p.y / Rh);\r
 107                        temp =Math.abs(this.lat0) - ProjConstants.HALF_PI;\r
 108                        if (Math.abs(temp) > ProjConstants.EPSLN) {\r
 109                           temp = cos_z -this.sin_lat_o * Math.sin(lat);\r
 110                           if(temp!=0.0) lon=ProjConstants.adjust_lon(this.long0+Math.atan2(p.x*sin_z*this.cos_lat_o,temp*Rh));\r
 111                        } else if (this.lat0 < 0.0) {\r
 112                           lon = ProjConstants.adjust_lon(this.long0 - Math.atan2(-p.x,p.y));\r
 113                        } else {\r
 114                           lon = ProjConstants.adjust_lon(this.long0 + Math.atan2(p.x, -p.y));\r
 115                        }\r
 116                     } else {\r
 117                       lat = this.lat0;\r
 118                     }\r
 119                     //return(OK);\r
 120                     p.x = lon;\r
 121                     p.y = lat;\r
 122                     return p;\r
 123                   }//lamazInv()\r
 124                 \r
 125                 \r
 126         }\r
 127 }