001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection;
003
004import java.util.Collections;
005import java.util.HashMap;
006import java.util.Map;
007import java.util.function.Consumer;
008import java.util.function.DoubleUnaryOperator;
009
010import org.openstreetmap.josm.data.Bounds;
011import org.openstreetmap.josm.data.ProjectionBounds;
012import org.openstreetmap.josm.data.coor.EastNorth;
013import org.openstreetmap.josm.data.coor.ILatLon;
014import org.openstreetmap.josm.data.coor.LatLon;
015import org.openstreetmap.josm.data.projection.datum.Datum;
016import org.openstreetmap.josm.data.projection.proj.Proj;
017import org.openstreetmap.josm.tools.Utils;
018
019/**
020 * Implementation of the Projection interface that represents a coordinate reference system and delegates
021 * the real projection and datum conversion to other classes.
022 *
023 * It handles false easting and northing, central meridian and general scale factor before calling the
024 * delegate projection.
025 *
026 * Forwards lat/lon values to the real projection in units of radians.
027 *
028 * The fields are named after Proj.4 parameters.
029 *
030 * Subclasses of AbstractProjection must set ellps and proj to a non-null value.
031 * In addition, either datum or nadgrid has to be initialized to some value.
032 */
033public abstract class AbstractProjection implements Projection {
034
035    protected Ellipsoid ellps;
036    protected Datum datum;
037    protected Proj proj;
038    protected double x0;            /* false easting (in meters) */
039    protected double y0;            /* false northing (in meters) */
040    protected double lon0;          /* central meridian */
041    protected double pm;            /* prime meridian */
042    protected double k0 = 1.0;      /* general scale factor */
043    protected double toMeter = 1.0; /* switch from meters to east/north coordinate units */
044
045    private volatile ProjectionBounds projectionBoundsBox;
046
047    /**
048     * Get the base ellipsoid that this projection uses.
049     * @return The {@link Ellipsoid}
050     */
051    public final Ellipsoid getEllipsoid() {
052        return ellps;
053    }
054
055    /**
056     * Gets the datum this projection is based on.
057     * @return The datum
058     */
059    public final Datum getDatum() {
060        return datum;
061    }
062
063    /**
064     * Replies the projection (in the narrow sense)
065     * @return The projection object
066     */
067    public final Proj getProj() {
068        return proj;
069    }
070
071    /**
072     * Gets an east offset that gets applied when converting the coordinate
073     * @return The offset to apply in meter
074     */
075    public final double getFalseEasting() {
076        return x0;
077    }
078
079    /**
080     * Gets an north offset that gets applied when converting the coordinate
081     * @return The offset to apply in meter
082     */
083    public final double getFalseNorthing() {
084        return y0;
085    }
086
087    /**
088     * Gets the meridian that this projection is centered on.
089     * @return The longitude of the meridian.
090     */
091    public final double getCentralMeridian() {
092        return lon0;
093    }
094
095    public final double getScaleFactor() {
096        return k0;
097    }
098
099    /**
100     * Get the factor that converts meters to intended units of east/north coordinates.
101     *
102     * For projected coordinate systems, the semi-major axis of the ellipsoid is
103     * always given in meters, which means the preliminary projection result will
104     * be in meters as well. This factor is used to convert to the intended units
105     * of east/north coordinates (e.g. feet in the US).
106     *
107     * For geographic coordinate systems, the preliminary "projection" result will
108     * be in degrees, so there is no reason to convert anything and this factor
109     * will by 1 by default.
110     *
111     * @return factor that converts meters to intended units of east/north coordinates
112     */
113    public final double getToMeter() {
114        return toMeter;
115    }
116
117    @Override
118    public EastNorth latlon2eastNorth(ILatLon toConvert) {
119        // TODO: Use ILatLon in datum, so we don't need to wrap it here.
120        LatLon ll = datum.fromWGS84(new LatLon(toConvert));
121        double[] en = proj.project(Utils.toRadians(ll.lat()), Utils.toRadians(LatLon.normalizeLon(ll.lon() - lon0 - pm)));
122        return new EastNorth(
123                (ellps.a * k0 * en[0] + x0) / toMeter,
124                (ellps.a * k0 * en[1] + y0) / toMeter);
125    }
126
127    @Override
128    public LatLon eastNorth2latlon(EastNorth en) {
129        // We know it is a latlon. Nice would be to change this method return type to ILatLon
130        return eastNorth2latlon(en, LatLon::normalizeLon);
131    }
132
133    @Override
134    public LatLon eastNorth2latlonClamped(EastNorth en) {
135        ILatLon ll = eastNorth2latlon(en, lon -> Utils.clamp(lon, -180, 180));
136        Bounds bounds = getWorldBoundsLatLon();
137        return new LatLon(Utils.clamp(ll.lat(), bounds.getMinLat(), bounds.getMaxLat()),
138                          Utils.clamp(ll.lon(), bounds.getMinLon(), bounds.getMaxLon()));
139    }
140
141    private LatLon eastNorth2latlon(EastNorth en, DoubleUnaryOperator normalizeLon) {
142        double[] latlonRad = proj.invproject(
143                 (en.east() * toMeter - x0) / ellps.a / k0,
144                (en.north() * toMeter - y0) / ellps.a / k0);
145        double lon = Utils.toDegrees(latlonRad[1]) + lon0 + pm;
146        LatLon ll = new LatLon(Utils.toDegrees(latlonRad[0]), normalizeLon.applyAsDouble(lon));
147        return datum.toWGS84(ll);
148    }
149
150    @Override
151    public Map<ProjectionBounds, Projecting> getProjectingsForArea(ProjectionBounds area) {
152        if (proj.lonIsLinearToEast()) {
153            //FIXME: Respect datum?
154            // wrap the world around
155            Bounds bounds = getWorldBoundsLatLon();
156            double minEast = latlon2eastNorth(bounds.getMin()).east();
157            double maxEast = latlon2eastNorth(bounds.getMax()).east();
158            double dEast = maxEast - minEast;
159            if ((area.minEast < minEast || area.maxEast > maxEast) && dEast > 0) {
160                // We could handle the dEast < 0 case but we don't need it atm.
161                int minChunk = (int) Math.floor((area.minEast - minEast) / dEast);
162                int maxChunk = (int) Math.floor((area.maxEast - minEast) / dEast);
163                HashMap<ProjectionBounds, Projecting> ret = new HashMap<>();
164                for (int chunk = minChunk; chunk <= maxChunk; chunk++) {
165                    ret.put(new ProjectionBounds(Math.max(area.minEast, minEast + chunk * dEast), area.minNorth,
166                                                 Math.min(area.maxEast, maxEast + chunk * dEast), area.maxNorth),
167                            new ShiftedProjecting(this, new EastNorth(-chunk * dEast, 0)));
168                }
169                return ret;
170            }
171        }
172
173        return Collections.singletonMap(area, this);
174    }
175
176    @Override
177    public double getDefaultZoomInPPD() {
178        // this will set the map scaler to about 1000 m
179        return 10 / getMetersPerUnit();
180    }
181
182    /**
183     * Returns The EPSG Code of this CRS.
184     * @return The EPSG Code of this CRS, null if it doesn't have one.
185     */
186    public abstract Integer getEpsgCode();
187
188    /**
189     * Default implementation of toCode().
190     * Should be overridden, if there is no EPSG code for this CRS.
191     */
192    @Override
193    public String toCode() {
194        return "EPSG:" + getEpsgCode();
195    }
196
197    protected static final double convertMinuteSecond(double minute, double second) {
198        return (minute/60.0) + (second/3600.0);
199    }
200
201    protected static final double convertDegreeMinuteSecond(double degree, double minute, double second) {
202        return degree + (minute/60.0) + (second/3600.0);
203    }
204
205    @Override
206    public final ProjectionBounds getWorldBoundsBoxEastNorth() {
207        ProjectionBounds result = projectionBoundsBox;
208        if (result == null) {
209            synchronized (this) {
210                result = projectionBoundsBox;
211                if (result == null) {
212                    ProjectionBounds bds = new ProjectionBounds();
213                    visitOutline(getWorldBoundsLatLon(), 1000, bds::extend);
214                    projectionBoundsBox = bds;
215                }
216            }
217        }
218        return projectionBoundsBox;
219    }
220
221    @Override
222    public Projection getBaseProjection() {
223        return this;
224    }
225
226    @Override
227    public void visitOutline(Bounds b, Consumer<EastNorth> visitor) {
228        visitOutline(b, 100, visitor);
229    }
230
231    private void visitOutline(Bounds b, int nPoints, Consumer<EastNorth> visitor) {
232        double maxlon = b.getMaxLon();
233        if (b.crosses180thMeridian()) {
234            maxlon += 360.0;
235        }
236        double spanLon = maxlon - b.getMinLon();
237        double spanLat = b.getMaxLat() - b.getMinLat();
238
239        //TODO: Use projection to see if there is any need for doing this along each axis.
240        for (int step = 0; step < nPoints; step++) {
241            visitor.accept(latlon2eastNorth(
242                    new LatLon(b.getMinLat(), b.getMinLon() + spanLon * step / nPoints)));
243        }
244        for (int step = 0; step < nPoints; step++) {
245            visitor.accept(latlon2eastNorth(
246                    new LatLon(b.getMinLat() + spanLat * step / nPoints, maxlon)));
247        }
248        for (int step = 0; step < nPoints; step++) {
249            visitor.accept(latlon2eastNorth(
250                    new LatLon(b.getMaxLat(), maxlon - spanLon * step / nPoints)));
251        }
252        for (int step = 0; step < nPoints; step++) {
253            visitor.accept(latlon2eastNorth(
254                    new LatLon(b.getMaxLat() - spanLat * step / nPoints, b.getMinLon())));
255        }
256    }
257}