Stereogram.java
001 /*
002
003 极射赤面投影(Stereogram projection)
004
005 PACKAGE: cma.common.projection
006 FILENAME: Coordinate.java
007 LANGUAGE: Java2 v1.4
008 ORIGINAL: none
009 DESCRIPTION:
010 CREATE: 2007-07-19 09:29:06
011 UPDATE: 2007-07-23
012 AUTHOR: 刘泽军 (BJ0773@gmail.com)
013 广西气象减灾研究所
014 Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
015 REFERENCE: http://mathworld.wolfram.com/StereographicProjection.html
016
017 Compile: javac Coordinate.java Stereogram.java
018
019 */
020
021 package cma.common.projection;
022
023 import java.io.*;
024 import java.awt.*;
025 import java.awt.geom.*;
026 import java.util.*;
027 import java.lang.Math.*;
028 import java.text.DecimalFormat;
029
030 import cma.common.algorithm.Geometry; //矢量的加法、减法、叉乘
031
032 public class Stereogram extends Coordinate {
033
034 /**
035 * 功能:
036 * 重置参数
037 * 参数:
038 * lon0,lat0 - 标准经纬度,lat0只能为90或-90,即BBQ或NBQ
039 * lon1,lat1 - 中心经纬度(用于定位,不一定是中心)
040 * x,y - 中心经纬度对应的屏幕坐标
041 * sc - 缩放系数
042 * 返回值:
043 * 无
044 */
045 private void reset ( //Degree
046 double lon0, double lat0,
047 double lon1, double lat1,
048 int x, int y,
049 double sc ) {
050 type = lat0 >= 0.0 ? Coordinate.BBQ : Coordinate.NBQ;
051 standard = new Point2D.Double ( Math.IEEEremainder ( Math.abs ( lon0 ) , 360.0 ) , lat0>= 0.0 ? 90.0 :- 90.0 ) ;
052 center = new Point2D.Double (
053 Math.IEEEremainder ( Math.abs ( lon1 ) , 360.0 ) ,
054 type == Coordinate.BBQ ?
055 ( lat1 < - 85.0 ? - 85.0 : lat1 > 90.0 ? 90.0 : lat1 ) : //北半球可以以北极为定位点
056 ( lat1 > 85.0 ? 85.0 : lat1 < - 90.0 ? - 90.0 : lat1 ) //南半球可以以南极为定位点
057 ) ;
058 place = new Point ( x, y ) ;
059 scale = sc <= 0.0 ? 1.0 : sc;
060 double _lon0 = Math.toRadians ( standard.x ) ;
061 double _lat0 = Math.toRadians ( standard.y ) ;
062 double _lon1 = Math.toRadians ( lon1 ) ;
063 double _lat1 = Math.toRadians ( lat1 ) ;
064 //0.04149是与Micaps1.0匹配的参数
065 double k = scale * 0.04149 * 2.0 *Coordinate.RADIUS/ ( 1.0 +Math.sin ( _lat0 ) *Math.sin ( _lat1 ) +Math.cos ( _lat0 ) *Math.cos ( _lat1 ) *Math.cos ( _lon1-_lon0 )) ;
066 double x0 = k * Math.cos ( _lat1 ) *Math.sin ( _lon1-_lon0 ) ;
067 double y0 = k * ( Math.cos ( _lat0 ) *Math.sin ( _lat1 ) -Math.sin ( _lat0 ) *Math.cos ( _lat1 ) *Math.cos ( _lon1-_lon0 )) ;
068 offset = new Point (( int )( 0.5 +place.x-x0 ) , ( int )( 0.5 +place.y-y0 )) ;
069 }
070
071 /**
072 * 功能:
073 * 构造函数
074 * 参数:
075 * 无(使用缺省值)
076 * 返回值:
077 * 无
078 */
079 public Stereogram () {
080 reset ( 110.0 , 90.0 , 109.40 , 24.35 , 640 , 480 , 1.0 ) ;
081 }
082
083 /**
084 * 功能:
085 * 构造函数
086 * 参数:
087 * lon,lat - 中心经纬度,
088 * px,py - 中心经纬度对应的屏幕坐标
089 * sc - 缩放系数
090 * 返回值:
091 * 无
092 */
093 public Stereogram ( double lon, double lat, int px, int py, double sc ) {
094 reset ( 110.0 , 35.0 , lon, lat, px, py, sc ) ;
095 }
096
097 /**
098 * 功能:
099 * 构造函数
100 * 参数:
101 * lon,lat - 中心经纬度,
102 * px,py - 中心经纬度对应的屏幕坐标
103 * sx,sy - 缩放比例
104 * 返回值:
105 * 无
106 */
107 public Stereogram ( double lon0, double lat0, double lon1, double lat1, int px, int py, double sc ) {
108 reset ( lon0, lat0, lon1, lat1, px, py, sc ) ;
109 }
110
111 /**
112 * 功能:
113 * 获得屏幕坐标
114 * 参数:
115 * lon - 经度
116 * lat - 纬度
117 * 返回值:
118 * 屏幕坐标
119 */
120 public Point getPosition ( double lon, double lat ) {
121 /*
122 double _lon0 = Math.toRadians(standard.x);
123 double _lat0 = Math.toRadians(standard.y);
124 double _lon1 = Math.toRadians(lon);
125 double _lat1 = Math.toRadians(lat);
126 //0.04149是与Micaps1.0匹配的参数,
127 double k = scale * 0.04149 * 2.0*Coordinate.RADIUS/(1.0+Math.sin(_lat0)*Math.sin(_lat1)+Math.cos(_lat0)*Math.cos(_lat1)*Math.cos(_lon1-_lon0));
128 double x = k * Math.cos(_lat1)*Math.sin(_lon1-_lon0);
129 double y = k * (Math.cos(_lat0)*Math.sin(_lat1)-Math.sin(_lat0)*Math.cos(_lat1)*Math.cos(_lon1-_lon0));
130 return(new Point((int)(0.5+offset.x+x), (int)(0.5+offset.y-y)));
131 */
132 double deltaLon = Math.toRadians ( lon-standard.x ) ;
133 double radianLat = Math.toRadians ( lat ) ;
134 //sign=sin(+/-90)=+/-1, cos(90)=cos(-90)=0.0,优化,去掉与标准纬度(即两个极点)有关的三角函数
135 double sgn = standard.y == 90.0 ? 1.0 : - 1.0 ;
136 //0.04149是与Micaps1.0匹配的参数,
137 double k = scale * 0.04149 * 2.0 *Coordinate.RADIUS/ ( 1.0 +sgn*Math.sin ( radianLat )) ;
138 double x = k * Math.cos ( radianLat ) *Math.sin ( deltaLon ) ;
139 double y = k * ( -sgn*Math.cos ( radianLat ) *Math.cos ( deltaLon )) ;
140 return ( new Point (( int )( 0.5 +offset.x+x ) , ( int )( 0.5 +offset.y-y ))) ;
141 }
142
143 /**
144 * 功能:
145 * 获得屏幕坐标对应的经纬度
146 * 参数:
147 * x - 屏幕水平坐标
148 * y - 屏幕垂直坐标
149 * 返回值:
150 * 对应的经纬度
151 */
152 public Point2D.Double getCoordinate ( int x, int y ) {
153 Point p0 = this .getPosition ( 0.0 , standard.y ) ;
154 if ( p0.x == x && p0.y == y ) {
155 return ( new Point2D.Double ( 0.0 , standard.y )) ;
156 }
157 Point p1 = this .getPosition ( 0.0 , 0.0 ) ;
158 if ( p1.x == x && p1.y == y ) {
159 return ( new Point2D.Double ( 0.0 , 0.0 )) ;
160 }
161 //求未知点与极点、(0.0E,0.0N)的夹角
162 double x1 = p0.x;
163 double y1 = p0.y;
164 double x2 = p1.x;
165 double y2 = p1.y;
166 double x3 = x;
167 double y3 = y;
168 //增加临时变量,但节约了计算量,由加减法17次,乘除法12次,调用函数5次,优化为加减法11次,乘除法8次,调用函数4次
169 double x3_x2 = x3-x2;
170 double y3_y2 = y3-y2;
171 double aa = x3_x2*x3_x2+y3_y2*y3_y2;
172 double x3_x1 = x3-x1;
173 double y3_y1 = y3-y1;
174 double bb = x3_x1*x3_x1+y3_y1*y3_y1;
175 double b = Math.sqrt ( bb ) ;
176 double x1_x2 = x1-x2;
177 double y1_y2 = y1-y2;
178 double cc = x1_x2*x1_x2+y1_y2*y1_y2;
179 double c = Math.sqrt ( cc ) ;
180 double angle = Math.acos (( bb+cc-aa ) / ( 2 *b*c )) ; //余弦定理,值在0-180之间
181 double pos = Geometry.pos ( p1.x, p1.y, p0.x, p0.y, x, y ) ; //判断点相对于0度经线的位置
182 if ( pos > 0.0 ) { //位于180-360度区域
183 angle = 2.0 * Math.PI - angle;
184 }
185 double lon = Math.toDegrees ( angle ) ;
186
187 double deltaLon = Math.toRadians ( lon-standard.x ) ;
188 double sgn = standard.y == 90.0 ? 1.0 : - 1.0 ;
189 double k = ( offset.y - y ) / ( -sgn * scale * 0.04149 * 2.0 * Coordinate.RADIUS * Math.cos ( deltaLon )) ;
190 double lat = Math.toDegrees ( Math.asin (( 1 -sgn*k*k ) / ( 1 +k*k ))) ; //由getPosition方法的公式反推得到
191
192 return ( new Point2D.Double ( lon, lat )) ;
193 }
194
195 /**
196 * 功能:获得相对于标准经度的偏角(逆时针方向)
197 * 参数:
198 * degreeLon - 经度值
199 * degreeLat - 纬度值
200 * 返回值:
201 * 偏角
202 */
203 public double getAngle ( double degreeLon, double degreeLat ) {
204 double agl = 360.0 +standard.x-degreeLon;
205 while ( agl>= 360.0 ) {
206 agl -= 360.0 ;
207 }
208 while ( agl<=- 360.0 ) {
209 agl += 360.0 ;
210 }
211 return ( agl ) ;
212 }
213
214 /**
215 * 功能:
216 * 画经线、纬线
217 * 参数:
218 * g - 图形设备
219 * f - 字体
220 * c - 画线颜色
221 * inc_lon - 经线间隔
222 * inc_lat - 纬线间隔
223 * 返回值:
224 * 无
225 */
226 public void drawGridLine ( Graphics2D g, Font f, Color c, int inc_lon, int inc_lat ) {
227
228 DecimalFormat df = new DecimalFormat ( "0.##########" ) ;
229 String lonString = df.format ( inc_lon ) ;
230 String latString = df.format ( inc_lat ) ;
231 DecimalFormat dfLon = new DecimalFormat ( "0." + ( "##########" .substring ( 0 , lonString.indexOf ( "." ) > 0 ?lonString.length () -lonString.indexOf ( "." ) : 1 ))) ;
232 DecimalFormat dfLat = new DecimalFormat ( "0." + ( "##########" .substring ( 0 , latString.indexOf ( "." ) > 0 ?latString.length () -latString.indexOf ( "." ) : 1 ))) ;
233
234 Color saveColor = g.getColor () ;
235 Font saveFont = g.getFont () ;
236 g.setColor ( c ) ;
237 g.setFont ( null ==f?f: new Font ( "Times New Roman" , Font.PLAIN, 12 )) ;
238 FontMetrics fm = g.getFontMetrics () ;
239 String text;
240 byte tmpByte [] ;
241 int bytesWidth, bytesHeight = fm.getHeight () ;
242 Point pos0, pos1, pos2;
243 pos0 = getPosition ( standard.x, type==Coordinate.BBQ? 90.0 :- 90.0 ) ;
244 double rotateAngle;
245 if ( inc_lon > 0 ) { //经线间隔为正值
246 for ( double lon= 0.0 ;lon< 360.0 ;lon=lon+inc_lon ) { //画经线
247 pos1 = getPosition ( lon, 0.0 ) ; //在赤道标经度值
248 g.translate ( pos1.x, pos1.y ) ;
249 rotateAngle = this .getAngle ( lon, 0.0 ) ;
250 g.rotate ( Math.toRadians ( rotateAngle )) ; //偏角
251 text = dfLon.format ( lon ) ;
252 tmpByte = text.getBytes () ;
253 bytesWidth = fm.bytesWidth ( tmpByte, 0 , tmpByte.length ) ;
254 g.drawString ( text, -bytesWidth/ 2 , bytesHeight/ 3 ) ;
255 g.rotate ( Math.toRadians ( -rotateAngle )) ; //偏角
256 g.translate ( -pos1.x, -pos1.y ) ;
257
258 pos1 = getPosition ( lon,- 90.0 + ( inc_lat> 0.0 ?inc_lat: 10.0 )) ;
259 pos2 = getPosition ( lon, 90.0 - ( inc_lat> 0.0 ?inc_lat: 10.0 )) ;
260 g.drawLine ( pos1.x, pos1.y, pos2.x, pos2.y ) ; //画经线
261 }
262 }
263 if ( inc_lat > 0 ) { //纬线间隔为正值
264 int r = 0 ;
265 for ( double lat= 90.0 -inc_lat;lat>=- 90.0 +inc_lat;lat=lat-inc_lat ) { //画纬线(纬圈)
266 pos1 = getPosition ( standard.x, lat ) ;
267 r = Math.abs ( pos1.y-pos0.y ) ;
268 g.drawArc ( pos0.x-r, pos0.y-r, 2 *r, 2 *r, 0 , 360 ) ;
269 }
270 for ( double lon= 0.0 ;lon< 360.0 ;lon=lon+ 90.0 ) { //在0、90、180、270经线上标纬度值
271 for ( double lat= 90.0 -inc_lat;lat>=- 90.0 +inc_lat;lat=lat-inc_lat ) {
272 pos1 = getPosition ( lon,lat ) ;
273 g.translate ( pos1.x, pos1.y ) ;
274 rotateAngle = this .getAngle ( lon, lat ) ;
275 g.rotate ( Math.toRadians ( rotateAngle )) ; //偏角
276 text = dfLat.format ( lat ) ;
277 tmpByte = text.getBytes () ;
278 bytesWidth = fm.bytesWidth ( tmpByte, 0 , tmpByte.length ) ;
279 g.drawString ( text, -bytesWidth/ 2 , bytesHeight/ 3 ) ;
280 // g.drawString(text, pos1.x-bytesWidth/2, pos1.y+bytesHeight/3);
281 g.rotate ( Math.toRadians ( -rotateAngle )) ; //偏角
282 g.translate ( -pos1.x, -pos1.y ) ;
283 }
284 }
285 }
286 g.setFont ( saveFont ) ;
287 g.setColor ( saveColor ) ;
288 }
289 }