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 }