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((ellps.a * k0 * en[0] + x0) / toMeter, (ellps.a * k0 * en[1] + y0) / toMeter);
123    }
124
125    @Override
126    public LatLon eastNorth2latlon(EastNorth en) {
127        // We know it is a latlon. Nice would be to change this method return type to ILatLon
128        return eastNorth2latlon(en, LatLon::normalizeLon);
129    }
130
131    @Override
132    public LatLon eastNorth2latlonClamped(EastNorth en) {
133        ILatLon ll = eastNorth2latlon(en, lon -> Utils.clamp(lon, -180, 180));
134        Bounds bounds = getWorldBoundsLatLon();
135        return new LatLon(Utils.clamp(ll.lat(), bounds.getMinLat(), bounds.getMaxLat()),
136                Utils.clamp(ll.lon(), bounds.getMinLon(), bounds.getMaxLon()));
137    }
138
139    private LatLon eastNorth2latlon(EastNorth en, DoubleUnaryOperator normalizeLon) {
140        double[] latlonRad = proj.invproject((en.east() * toMeter - x0) / ellps.a / k0, (en.north() * toMeter - y0) / ellps.a / k0);
141        double lon = Utils.toDegrees(latlonRad[1]) + lon0 + pm;
142        LatLon ll = new LatLon(Utils.toDegrees(latlonRad[0]), normalizeLon.applyAsDouble(lon));
143        return datum.toWGS84(ll);
144    }
145
146    @Override
147    public Map<ProjectionBounds, Projecting> getProjectingsForArea(ProjectionBounds area) {
148        if (proj.lonIsLinearToEast()) {
149            //FIXME: Respect datum?
150            // wrap the world around
151            Bounds bounds = getWorldBoundsLatLon();
152            double minEast = latlon2eastNorth(bounds.getMin()).east();
153            double maxEast = latlon2eastNorth(bounds.getMax()).east();
154            double dEast = maxEast - minEast;
155            if ((area.minEast < minEast || area.maxEast > maxEast) && dEast > 0) {
156                // We could handle the dEast < 0 case but we don't need it atm.
157                int minChunk = (int) Math.floor((area.minEast - minEast) / dEast);
158                int maxChunk = (int) Math.floor((area.maxEast - minEast) / dEast);
159                HashMap<ProjectionBounds, Projecting> ret = new HashMap<>();
160                for (int chunk = minChunk; chunk <= maxChunk; chunk++) {
161                    ret.put(new ProjectionBounds(Math.max(area.minEast, minEast + chunk * dEast), area.minNorth,
162                            Math.min(area.maxEast, maxEast + chunk * dEast), area.maxNorth),
163                            new ShiftedProjecting(this, new EastNorth(-chunk * dEast, 0)));
164                }
165                return ret;
166            }
167        }
168
169        return Collections.singletonMap(area, this);
170    }
171
172    @Override
173    public double getDefaultZoomInPPD() {
174        // this will set the map scaler to about 1000 m
175        return 10 / getMetersPerUnit();
176    }
177
178    /**
179     * @return The EPSG Code of this CRS, null if it doesn't have one.
180     */
181    public abstract Integer getEpsgCode();
182
183    /**
184     * Default implementation of toCode().
185     * Should be overridden, if there is no EPSG code for this CRS.
186     */
187    @Override
188    public String toCode() {
189        return "EPSG:" + getEpsgCode();
190    }
191
192    protected static final double convertMinuteSecond(double minute, double second) {
193        return (minute/60.0) + (second/3600.0);
194    }
195
196    protected static final double convertDegreeMinuteSecond(double degree, double minute, double second) {
197        return degree + (minute/60.0) + (second/3600.0);
198    }
199
200    @Override
201    public final ProjectionBounds getWorldBoundsBoxEastNorth() {
202        ProjectionBounds result = projectionBoundsBox;
203        if (result == null) {
204            synchronized (this) {
205                result = projectionBoundsBox;
206                if (result == null) {
207                    ProjectionBounds bds = new ProjectionBounds();
208                    visitOutline(getWorldBoundsLatLon(), 1000, bds::extend);
209                    projectionBoundsBox = bds;
210                }
211            }
212        }
213        return projectionBoundsBox;
214    }
215
216    @Override
217    public Projection getBaseProjection() {
218        return this;
219    }
220
221    @Override
222    public void visitOutline(Bounds b, Consumer<EastNorth> visitor) {
223        visitOutline(b, 100, visitor);
224    }
225
226    private void visitOutline(Bounds b, int nPoints, Consumer<EastNorth> visitor) {
227        double maxlon = b.getMaxLon();
228        if (b.crosses180thMeridian()) {
229            maxlon += 360.0;
230        }
231        double spanLon = maxlon - b.getMinLon();
232        double spanLat = b.getMaxLat() - b.getMinLat();
233
234        //TODO: Use projection to see if there is any need for doing this along each axis.
235        for (int step = 0; step < nPoints; step++) {
236            visitor.accept(latlon2eastNorth(
237                    new LatLon(b.getMinLat(), b.getMinLon() + spanLon * step / nPoints)));
238        }
239        for (int step = 0; step < nPoints; step++) {
240            visitor.accept(latlon2eastNorth(
241                    new LatLon(b.getMinLat() + spanLat * step / nPoints, maxlon)));
242        }
243        for (int step = 0; step < nPoints; step++) {
244            visitor.accept(latlon2eastNorth(
245                    new LatLon(b.getMaxLat(), maxlon - spanLon * step / nPoints)));
246        }
247        for (int step = 0; step < nPoints; step++) {
248            visitor.accept(latlon2eastNorth(
249                    new LatLon(b.getMaxLat() - spanLat * step / nPoints, b.getMinLon())));
250        }
251    }
252}