Lambert(兰勃托)投影--我国天气图底图广泛采用的一种投影

---

Lambert.java

---

001  /**
002 
003     Lambert兰勃特投影
004 
005        PACKAGE: cma.common.projection
006       FILENAME: Lambert.java
007       LANGUAGE: Java2 v1.4
008       ORIGINAL: 最初由成气院向卫国袁东升老师编写，但与Micaps1.0的兰勃特投影不一致，距离标准经度远的地方误差增大。
009    DESCRIPTION: Lambert projection coordinate
010         CREATE: 2003-03-06 00:31:21
011         UPDATE: 2007-07-17 根据 http://mathworld.wolfram.com  提供的公式重新编写
012         AUTHOR: 刘泽军 (BJ0773@gmail.com)
013                 广西气象减灾研究所
014                 Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
015      REFERENCE: http://mathworld.wolfram.com/LambertConformalConicProjection.html
016 
017        COMPILE: javac Coordinate.java Lambert.java
018 
019  ********** 重写后，与Micaps1.0版的“T-兰勃脱投影”完全一致！ **********
020 
021        HISTORY:
022           0、lbt.cpp，两位老师的源代码
023           1、CCoordinate.cpp，线性和Lambert投影，主要用于FY-2B卫星云图定位和Micaps交互；
024           2、CRotateCoordinate.cpp，坐标旋转，主要是用于柳州Doppler雷达选址时从扫描的军用地图中定位及读取高程信息)；
025           3、CLambert.cpp、CLines.cpp，各投影代码独立为类
026           4、Lambert.java，兰勃托投影改写为JAVA版
027           5、Lambert.java，2007-07-17重新编写
028 
029  **/
030 
031  package  cma.common.projection;
032 
033  import  java.io.*;
034  import  java.awt.*;
035  import  java.awt.geom.*;
036  import  java.util.*;
037  import  java.lang.Math.*;
038  import  java.text.DecimalFormat;
039 
040  import  cma.common.atmos.*;
041 
042  public class  Lambert  extends  Coordinate  {
043 
044       //private   double  fn = 0.715566799999999;
045       private  double   F;
046       private  double   n;
047       private  double   p0;
048 
049       private  void  reset ( //Degree
050               double   standardLon,     double   standardLat,
051               double   positionLon,     double   positionLat,
052               int      positionHori,    int      positionVert,
053               double   zoomIndex ) {
054           type            = Coordinate.LAMBERT;
055           standard        =  new  Point2D.Double ( Math.IEEEremainder ( Math.abs ( standardLon ) ,  360.0 ) ,  0.0 ) ;
056           center          =  new  Point2D.Double ( Math.IEEEremainder ( Math.abs ( positionLon ) ,  360.0 ) , positionLat < - 85.0  ? - 85.0  : positionLat >  85.0  ?  85.0  : positionLat ) ;
057           place           =  new  Point ( positionHori, positionVert ) ;
058           scale           = zoomIndex <=  0.0  ?  1.0  : zoomIndex;
059           scaleOriginal   = scale;
060 
061           double   lat1= Math.toRadians ( 30.0 ) , lat2 = Math.toRadians ( 60.0 ) ; //标准纬度
062           double   lat3= Math.toRadians ( 45.0 + 30.0 / 2.0 ) , lat4 = Math.toRadians ( 45.0 + 60.0 / 2 ) ; //标准纬度
063           double   lat5= Math.toRadians ( 45.0 +center.y/ 2.0 ) ,  lat6 = Math.toRadians ( 45.0 +standard.y/ 2.0 ) ;
064           n   = Math.log ( Math.cos ( lat1 ) /Math.cos ( lat2 )) /Math.log ( Math.tan ( lat4 ) /Math.tan ( lat3 )) ;
065           F   = scale *  484.96  * Math.cos ( lat1 ) *Math.pow ( Math.tan ( lat3 ) ,n ) /n; //484.96为根据Micaps1.0拟合得到的数值
066           p0  = F * Math.pow ( 1.0 /Math.tan ( lat6 ) ,n ) ;
067           double   r   = Math.toRadians ( n* ( center.x-standard.x )) ;
068           double   p   = F * Math.pow ( 1.0 /Math.tan ( lat5 ) ,n ) ;
069           double   x   = p * Math.sin ( r ) ;
070           double   y   = p * Math.cos ( r )  - p0;
071           offset      =  new  Point (( int )( 0.5  + place.x - x ) ,  ( int )( 0.5  + place.y - y )) ;
072       }
073 
074       public  Lambert () {
075           reset (
076               110.0 ,   30.0 ,
077               110.0 ,   35.0 ,
078               512 ,     384 ,
079               1
080           ) ;
081       }
082 
083       public  Lambert ( //Degree
084               double   standardLon,     double   standardLat,
085               double   positionLon,     double   positionLat,
086               int      positionHori,    int      positionVert,
087               double   zoomIndex ) {
088           reset (
089               standardLon,standardLat,
090               positionLon,    positionLat,
091               positionHori,   positionVert,
092               zoomIndex
093           ) ;
094       }
095 
096       public  Lambert ( //Degree
097               double   positionLon,     double   positionLat,
098               int      positionHori,    int      positionVert,
099               double   zoomIndex ) {
100           reset (
101               110.0 ,       30.0 ,
102               positionLon,    positionLat,
103               positionHori,   positionVert,
104               zoomIndex
105           ) ;
106       }
107 
108  /**
109    * 功能：获得屏幕坐标
110    * 参数：
111    *      lon     - 经度值(度)
112    *      lat     - 纬度值(度)
113    * 返回值：
114    *      屏幕坐标
115    */
116       public  Point getPosition ( double  lon,  double  lat ) {
117           double   lat5= Math.toRadians ( 45.0 +lat/ 2.0 ) ;
118           double   r   = Math.toRadians ( n* ( lon-standard.x )) ;
119           double   p   = F * Math.pow ( 1.0 /Math.tan ( lat5 ) ,n ) ;
120           double   x   = p * Math.sin ( r ) ;
121           double   y   = p * Math.cos ( r )  - p0;
122           return ( new  Point (( int )( 0.5  + offset.x + x ) ,  ( int )( 0.5  + offset.y + y ))) ;
123       }
124 
125  /**
126    * 功能：获得经纬度坐标
127    * 参数：
128    *  xScreen     - 屏幕坐标值(x方向)
129    *  yScreen     - 屏幕坐标值(y方向)
130    * 返回值：
131    *      经纬度坐标
132    */
133       public  Point2D.Double getCoordinate ( int  x,  int  y ) {
134 
135           double   lon, lat, r;
136           if (  y-offset.y+p0 ==  0  ) {
137               lon = standard.x+ 45.0 /n;
138               lat =  2.0 * ( Math.toDegrees ( Math.atan ( 1.0 /Math.pow (( x-offset.x ) /F, 1.0 /n ))) - 45.0 ) ;
139           }
140           else  {
141               lon = standard.x + Math.toDegrees ( Math.atan (( x-offset.x ) / ( y-offset.y+p0 ))) /n;
142               r   = Math.toRadians ( n* ( lon-standard.x )) ;
143               lat =  2.0 * ( Math.toDegrees ( Math.atan ( 1.0 /Math.pow (( y-offset.y+p0 ) /Math.cos ( r ) /F, 1.0 /n ))) - 45.0 ) ;
144           }
145           return ( new  Point2D.Double ( lon,lat )) ;
146       }
147 
148  /**
149    * 功能：获得相对于标准经度的偏角
150    * 参数：
151    *  degreeLon   - 经度值
152    *  degreeLat   - 纬度值
153    * 返回值：
154    *      偏角
155    */
156       public  double  getAngle ( double  degreeLon,  double  degreeLat ) {
157           if (  9999.0  <= degreeLat || - 9999.0  >= degreeLat ||  90.0  == degreeLat || - 90.0  == degreeLat  ) {
158               return ( 0.0 ) ;
159           }
160           Point   xy1     = getPosition ( standard.x,  90.0 ) ;
161           Point   xy2     = getPosition ( degreeLon, degreeLat ) ;
162           double   x1      = xy1.x;
163           double   y1      = xy1.y;
164           double   x2      = xy2.x;
165           double   y2      = xy2.y;
166           double   rotateAngle =  0.0 ;
167           if (  x2 == x1  ) {
168               rotateAngle = y2 >= y1 ?  0.0  :  180.0 ;
169           }
170           else if (  y1 == y2  ) {
171               rotateAngle = x2 >= x1 ?  270.0  :  90.0 ;
172           }
173           else  { //分象限判断旋转角度，坐标原点在xy1
174               double   deltaY      = Math.abs ( y2 - y1 ) ;
175               double   r           = Math.sqrt (( x2 - x1 ) * ( x2 - x1 )  +  ( y2 - y1 ) * ( y2 - y1 )) ;
176               double   angelDegree = Math.toDegrees ( Math.asin ( deltaY/r )) ;
177               if (  x2 > x1 && y2 > y1  ) {
178                   rotateAngle = - 90.0  + angelDegree;
179               }
180               else if (  x2 > x1 && y2 < y1  ) {
181                   rotateAngle = - 90.0  - angelDegree;
182               }
183               else if (  x2 < x1 && y2 < y1  ) {
184                   rotateAngle =  90.0  + angelDegree;
185               }
186               else if (  x2 < x1 && y2 > y1  ) {
187                   rotateAngle =  90.0  - angelDegree;
188               }
189               else  {
190                   return ( Algorithm.DefaultValue ) ;
191               }
192           }
193           return ( rotateAngle ) ;
194       }
195 
196  /**
197    * 功能：
198    *      画经线、纬线
199    * 参数：
200    *      g       - 图形设备
201    *      f       - 字体
202    *      c       - 画线颜色
203    *      inc_lon - 经线间隔
204    *      inc_lat - 纬线间隔
205    * 返回值：
206    *      无
207    */
208       public  void  drawGridLine ( Graphics2D g, Font f, Color c,  int  inc_lon,  int  inc_lat ) {
209 
210           DecimalFormat   df  =  new  DecimalFormat ( "0.#" ) ;
211           Color   saveColor   = g.getColor () ;
212           Font    saveFont    = g.getFont () ;
213           g.setColor ( c ) ;
214           g.setFont ( null ==f?f: new  Font ( "Times New Roman" , Font.PLAIN,  12 )) ;
215           Point       xy, xy1;
216           double       rotateAngle;
217           FontMetrics fm  = g.getFontMetrics () ;
218           String      text;
219           byte         tmpByte [] ;
220           int          bytesWidth, bytesHeight = fm.getHeight () ;;
221           g.setColor ( c ) ;
222           Point   pos1, pos2;
223           //画纬线
224           for ( int  lat= 80 ;lat>=- 80 +inc_lat;lat=lat-inc_lat ) {
225               for ( int  lon= 0 ;lon< 360 ;lon++ ) {
226                   pos1    = getPosition ( lon, lat ) ;
227                   pos2    = getPosition ( lon+ 1 , lat ) ;
228                   g.drawLine ( pos1.x, pos1.y, pos2.x, pos2.y ) ;
229                   if (  lon %  90  ==  0  ) {
230                       g.translate ( pos1.x, pos1.y ) ;
231                       rotateAngle = getAngle ( lon, lat ) ;
232                       if (  rotateAngle != Algorithm.DefaultValue  ) {
233                           g.rotate ( Math.toRadians ( rotateAngle )) ;   //偏角
234                       }
235                       text        = df.format ( lat ) ;
236                       tmpByte     = text.getBytes () ;
237                       bytesWidth  = fm.bytesWidth ( tmpByte,  0 , tmpByte.length ) ;
238                       g.drawString ( text, -bytesWidth/ 2 , bytesHeight/ 3 ) ; //标纬度值
239                       if (  rotateAngle != Algorithm.DefaultValue  ) {
240                           g.rotate ( Math.toRadians ( -rotateAngle )) ;  //偏角
241                       }
242                       g.translate ( -pos1.x, -pos1.y ) ;
243                   }
244               }
245           }
246           //画经线
247           for ( int  lon= 0 ;lon<= 360 ;lon=lon+inc_lon ) {
248               pos1    = getPosition ( lon,   80.0 ) ;
249               pos2    = getPosition ( lon, - 80.0 ) ;
250               g.drawLine ( pos1.x, pos1.y, pos2.x, pos2.y ) ;
251               pos1    = getPosition ( lon,    0.0 ) ;
252               g.translate ( pos1.x, pos1.y ) ;
253               rotateAngle = getAngle ( lon,  0.0 ) ;
254               if (  rotateAngle != Algorithm.DefaultValue  ) {
255                   g.rotate ( Math.toRadians ( rotateAngle )) ;   //偏角
256               }
257               text        = df.format ( lon ) ;
258               tmpByte     = text.getBytes () ;
259               bytesWidth  = fm.bytesWidth ( tmpByte,  0 , tmpByte.length ) ;
260               g.drawString ( text, -bytesWidth/ 2 , bytesHeight/ 3 ) ; //标纬度值
261               if (  rotateAngle != Algorithm.DefaultValue  ) {
262                   g.rotate ( Math.toRadians ( -rotateAngle )) ;  //偏角
263               }
264               g.translate ( -pos1.x, -pos1.y ) ;
265           }
266           g.setFont ( saveFont ) ;
267           g.setColor ( saveColor ) ;
268       }
269  }