001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.3_testing/src/org/deegree/framework/util/MapUtils.java $
002    /*----------------------------------------------------------------------------
003     This file is part of deegree, http://deegree.org/
004     Copyright (C) 2001-2009 by:
005     Department of Geography, University of Bonn
006     and
007     lat/lon GmbH
008    
009     This library is free software; you can redistribute it and/or modify it under
010     the terms of the GNU Lesser General Public License as published by the Free
011     Software Foundation; either version 2.1 of the License, or (at your option)
012     any later version.
013     This library is distributed in the hope that it will be useful, but WITHOUT
014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
016     details.
017     You should have received a copy of the GNU Lesser General Public License
018     along with this library; if not, write to the Free Software Foundation, Inc.,
019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
020    
021     Contact information:
022    
023     lat/lon GmbH
024     Aennchenstr. 19, 53177 Bonn
025     Germany
026     http://lat-lon.de/
027    
028     Department of Geography, University of Bonn
029     Prof. Dr. Klaus Greve
030     Postfach 1147, 53001 Bonn
031     Germany
032     http://www.geographie.uni-bonn.de/deegree/
033    
034     e-mail: info@deegree.org
035     ----------------------------------------------------------------------------*/
036    package org.deegree.framework.util;
037    
038    import static java.lang.Math.sqrt;
039    
040    import org.deegree.framework.log.ILogger;
041    import org.deegree.framework.log.LoggerFactory;
042    import org.deegree.i18n.Messages;
043    import org.deegree.model.crs.CRSFactory;
044    import org.deegree.model.crs.CoordinateSystem;
045    import org.deegree.model.crs.GeoTransformer;
046    import org.deegree.model.spatialschema.Envelope;
047    import org.deegree.model.spatialschema.GeometryFactory;
048    import org.deegree.model.spatialschema.Position;
049    
050    /**
051     * 
052     * 
053     * @version $Revision: 20692 $
054     * @author <a href="mailto:poth@lat-lon.de">Andreas Poth</a>
055     * @author last edited by: $Author: apoth $
056     * 
057     * @version 1.0. $Revision: 20692 $, $Date: 2009-11-10 13:49:55 +0100 (Di, 10. Nov 2009) $
058     * 
059     * @since 2.0
060     */
061    public class MapUtils {
062    
063        private static ILogger LOG = LoggerFactory.getLogger( MapUtils.class );
064    
065        /**
066         * The value of sqrt(2)
067         */
068        public static final double SQRT2 = sqrt( 2 );
069    
070        /**
071         * The Value of a PixelSize
072         */
073        public static final double DEFAULT_PIXEL_SIZE = 0.00028;
074    
075        /**
076         * @param mapWidth
077         * @param mapHeight
078         * @param bbox
079         * @param crs
080         * @return the WMS 1.1.1 scale (size of the diagonal pixel)
081         */
082        public static double calcScaleWMS111( int mapWidth, int mapHeight, Envelope bbox, CoordinateSystem crs ) {
083            if ( mapWidth == 0 || mapHeight == 0 ) {
084                return 0;
085            }
086            double scale = 0;
087    
088            if ( crs == null ) {
089                throw new RuntimeException( "Invalid crs: " + crs );
090            }
091    
092            try {
093                if ( "m".equalsIgnoreCase( crs.getAxisUnits()[0].toString() ) ) {
094                    /*
095                     * this method to calculate a maps scale as defined in OGC WMS and SLD specification is not required for
096                     * maps having a projected reference system. Direct calculation of scale avoids uncertainties
097                     */
098                    double dx = bbox.getWidth() / mapWidth;
099                    double dy = bbox.getHeight() / mapHeight;
100                    scale = sqrt( dx * dx + dy * dy );
101                } else {
102    
103                    if ( !crs.getIdentifier().equalsIgnoreCase( "EPSG:4326" ) ) {
104                        // transform the bounding box of the request to EPSG:4326
105                        GeoTransformer trans = new GeoTransformer( CRSFactory.create( "EPSG:4326" ) );
106                        bbox = trans.transform( bbox, crs );
107                    }
108                    double dx = bbox.getWidth() / mapWidth;
109                    double dy = bbox.getHeight() / mapHeight;
110                    Position min = GeometryFactory.createPosition( bbox.getMin().getX() + dx * ( mapWidth / 2d - 1 ),
111                                                                   bbox.getMin().getY() + dy * ( mapHeight / 2d - 1 ) );
112                    Position max = GeometryFactory.createPosition( bbox.getMin().getX() + dx * ( mapWidth / 2d ),
113                                                                   bbox.getMin().getY() + dy * ( mapHeight / 2d ) );
114    
115                    double distance = calcDistance( min.getX(), min.getY(), max.getX(), max.getY() );
116    
117                    scale = distance / SQRT2;
118    
119                }
120            } catch ( Exception e ) {
121                LOG.logError( e.getMessage(), e );
122                throw new RuntimeException( Messages.getMessage( "FRAMEWORK_ERROR_SCALE_CALC", e.getMessage() ) );
123            }
124    
125            return scale;
126        }
127    
128        /**
129         * @param mapWidth
130         * @param mapHeight
131         * @param bbox
132         * @param crs
133         * @return the WMS 1.3.0 scale (horizontal size of the pixel, pixel size == 0.28mm)
134         */
135        public static double calcScaleWMS130( int mapWidth, int mapHeight, Envelope bbox, CoordinateSystem crs ) {
136            if ( mapWidth == 0 || mapHeight == 0 ) {
137                return 0;
138            }
139    
140            double scale = 0;
141    
142            if ( crs == null ) {
143                throw new RuntimeException( "Invalid crs: " + crs );
144            }
145    
146            try {
147                if ( "m".equalsIgnoreCase( crs.getAxisUnits()[0].toString() ) ) {
148                    /*
149                     * this method to calculate a maps scale as defined in OGC WMS and SLD specification is not required for
150                     * maps having a projected reference system. Direct calculation of scale avoids uncertainties
151                     */
152                    double dx = bbox.getWidth() / mapWidth;
153                    scale = dx / DEFAULT_PIXEL_SIZE;
154                } else {
155    
156                    if ( !crs.getIdentifier().equalsIgnoreCase( "EPSG:4326" ) ) {
157                        // transform the bounding box of the request to EPSG:4326
158                        GeoTransformer trans = new GeoTransformer( CRSFactory.create( "EPSG:4326" ) );
159                        bbox = trans.transform( bbox, crs );
160                    }
161                    double dx = bbox.getWidth() / mapWidth;
162                    double dy = bbox.getHeight() / mapHeight;
163    
164                    double minx = bbox.getMin().getX() + dx * ( mapWidth / 2d - 1 );
165                    double miny = bbox.getMin().getY() + dy * ( mapHeight / 2d - 1 );
166                    double maxx = bbox.getMin().getX() + dx * ( mapWidth / 2d );
167                    double maxy = bbox.getMin().getY() + dy * ( mapHeight / 2d - 1 );
168    
169                    double distance = calcDistance( minx, miny, maxx, maxy );
170    
171                    scale = distance / SQRT2 / DEFAULT_PIXEL_SIZE;
172    
173                }
174            } catch ( Exception e ) {
175                LOG.logError( e.getMessage(), e );
176                throw new RuntimeException( Messages.getMessage( "FRAMEWORK_ERROR_SCALE_CALC", e.getMessage() ) );
177            }
178    
179            return scale;
180        }
181    
182        /**
183         * calculates the map scale (denominator) as defined in the OGC SLD 1.0.0 specification
184         * 
185         * @param mapWidth
186         *            map width in pixel
187         * @param mapHeight
188         *            map height in pixel
189         * @param bbox
190         *            bounding box of the map
191         * @param crs
192         *            coordinate reference system of the map
193         * @param pixelSize
194         *            size of one pixel of the map measured in meter
195         * 
196         * @return a maps scale based on the diagonal size of a pixel at the center of the map in meter.
197         * @throws RuntimeException
198         */
199        public static double calcScale( int mapWidth, int mapHeight, Envelope bbox, CoordinateSystem crs, double pixelSize )
200                                throws RuntimeException {
201    
202            double sqpxsize;
203            if ( pixelSize == 1d ) {
204                LOG.logDebug( "Calculating WMS 1.1.1 scale." );
205                return calcScaleWMS111( mapWidth, mapHeight, bbox, crs );
206            } else if ( pixelSize == DEFAULT_PIXEL_SIZE ) {
207                LOG.logDebug( "Calculating WMS 1.3.0 scale." );
208                return calcScaleWMS130( mapWidth, mapHeight, bbox, crs );
209            } else {
210                sqpxsize = pixelSize * pixelSize;
211                sqpxsize += sqpxsize;
212                sqpxsize = sqrt( sqpxsize );
213            }
214    
215            if ( mapWidth == 0 || mapHeight == 0 ) {
216                return 0;
217            }
218    
219            double scale = 0;
220    
221            CoordinateSystem cs = crs;
222    
223            if ( cs == null ) {
224                throw new RuntimeException( "Invalid crs: " + crs );
225            }
226    
227            try {
228                if ( "m".equalsIgnoreCase( cs.getAxisUnits()[0].toString() ) ) {
229                    /*
230                     * this method to calculate a maps scale as defined in OGC WMS and SLD specification is not required for
231                     * maps having a projected reference system. Direct calculation of scale avoids uncertainties
232                     */
233                    double dx = bbox.getWidth() / mapWidth;
234                    double dy = bbox.getHeight() / mapHeight;
235                    scale = Math.sqrt( dx * dx + dy * dy ) / sqpxsize;
236                } else {
237    
238                    if ( !crs.getIdentifier().equalsIgnoreCase( "EPSG:4326" ) ) {
239                        // transform the bounding box of the request to EPSG:4326
240                        GeoTransformer trans = new GeoTransformer( CRSFactory.create( "EPSG:4326" ) );
241                        bbox = trans.transform( bbox, crs );
242                    }
243                    double dx = bbox.getWidth() / mapWidth;
244                    double dy = bbox.getHeight() / mapHeight;
245                    Position min = GeometryFactory.createPosition( bbox.getMin().getX() + dx * ( mapWidth / 2d - 1 ),
246                                                                   bbox.getMin().getY() + dy * ( mapHeight / 2d - 1 ) );
247                    Position max = GeometryFactory.createPosition( bbox.getMin().getX() + dx * ( mapWidth / 2d ),
248                                                                   bbox.getMin().getY() + dy * ( mapHeight / 2d ) );
249    
250                    double distance = calcDistance( min.getX(), min.getY(), max.getX(), max.getY() );
251    
252                    scale = distance / sqpxsize;
253    
254                }
255            } catch ( Exception e ) {
256                LOG.logError( e.getMessage(), e );
257                throw new RuntimeException( Messages.getMessage( "FRAMEWORK_ERROR_SCALE_CALC", e.getMessage() ) );
258            }
259    
260            return scale;
261    
262        }
263    
264        /**
265         * calculates the distance in meters between two points in EPSG:4326 coodinates. this is a convenience method
266         * assuming the world is a ball
267         * 
268         * @param lon1
269         * @param lat1
270         * @param lon2
271         * @param lat2
272         * @return the distance in meters between two points in EPSG:4326 coords
273         */
274        public static double calcDistance( double lon1, double lat1, double lon2, double lat2 ) {
275            double r = 6378.137;
276            double rad = Math.PI / 180d;
277            double cose = Math.sin( rad * lon1 ) * Math.sin( rad * lon2 ) + Math.cos( rad * lon1 ) * Math.cos( rad * lon2 )
278                          * Math.cos( rad * ( lat1 - lat2 ) );
279            double dist = r * Math.acos( cose ) * Math.cos( rad * Math.min( lat1, lat2 ) );
280    
281            // * 0.835 is just an heuristic correction factor
282            return dist * 1000 * 0.835;
283        }
284    
285        /**
286         * The method calculates a new Envelope from the <code>requestedBarValue</code> It will either zoom in or zoom out
287         * of the <code>actualBBOX<code> depending
288         * on the ratio of the <code>requestedBarValue</code> to the <code>actualBarValue</code>
289         * 
290         * @param currentEnvelope
291         *            current Envelope
292         * @param currentScale
293         *            the scale of the current envelope
294         * @param requestedScale
295         *            requested scale value
296         * @return a new Envelope
297         */
298        public static Envelope scaleEnvelope( Envelope currentEnvelope, double currentScale, double requestedScale ) {
299    
300            double ratio = requestedScale / currentScale;
301            double newWidth = currentEnvelope.getWidth() * ratio;
302            double newHeight = currentEnvelope.getHeight() * ratio;
303            double midX = currentEnvelope.getMin().getX() + ( currentEnvelope.getWidth() / 2d );
304            double midY = currentEnvelope.getMin().getY() + ( currentEnvelope.getHeight() / 2d );
305    
306            double minx = midX - newWidth / 2d;
307            double maxx = midX + newWidth / 2d;
308            double miny = midY - newHeight / 2d;
309            double maxy = midY + newHeight / 2d;
310    
311            return GeometryFactory.createEnvelope( minx, miny, maxx, maxy, currentEnvelope.getCoordinateSystem() );
312    
313        }
314    
315        /**
316         * This method ensures the bbox is resized (shrunk) to match the aspect ratio defined by mapHeight/mapWidth
317         * 
318         * @param bbox
319         * @param mapWith
320         * @param mapHeight
321         * @return a new bounding box with the aspect ratio given my mapHeight/mapWidth
322         */
323        public static final Envelope ensureAspectRatio( Envelope bbox, double mapWith, double mapHeight ) {
324    
325            double minx = bbox.getMin().getX();
326            double miny = bbox.getMin().getY();
327            double maxx = bbox.getMax().getX();
328            double maxy = bbox.getMax().getY();
329    
330            double dx = maxx - minx;
331            double dy = maxy - miny;
332    
333            double ratio = mapHeight / mapWith;
334    
335            if ( dx >= dy ) {
336                // height has to be corrected
337                double[] normCoords = getNormalizedCoords( dx, ratio, miny, maxy );
338                miny = normCoords[0];
339                maxy = normCoords[1];
340            } else {
341                // width has to be corrected
342                ratio = mapWith / mapHeight;
343                double[] normCoords = getNormalizedCoords( dy, ratio, minx, maxx );
344                minx = normCoords[0];
345                maxx = normCoords[1];
346            }
347            CoordinateSystem crs = bbox.getCoordinateSystem();
348    
349            return GeometryFactory.createEnvelope( minx, miny, maxx, maxy, crs );
350        }
351    
352        private static final double[] getNormalizedCoords( double normLen, double ratio, double min, double max ) {
353            double mid = ( max - min ) / 2 + min;
354            min = mid - ( normLen / 2 ) * ratio;
355            max = mid + ( normLen / 2 ) * ratio;
356            double[] newCoords = { min, max };
357            return newCoords;
358        }
359    
360    }