Project

General

Profile

Download (8.46 KB) Statistics
| Branch: | Tag: | Revision:
1
/**
2
 * Copyright (C) 2020 EDIT
3
 * European Distributed Institute of Taxonomy
4
 * http://www.e-taxonomy.eu
5
 *
6
 * The contents of this file are subject to the Mozilla Public License Version 1.1
7
 * See LICENSE.TXT at the top of this package for the full license terms.
8
 */
9

    
10
package eu.etaxonomy.cdm.ext.geo.kml;
11

    
12
import java.awt.geom.Point2D;
13

    
14
import javax.measure.Quantity;
15
import javax.measure.Unit;
16
import javax.measure.UnitConverter;
17
import javax.measure.quantity.Length;
18

    
19
import org.apache.logging.log4j.LogManager;
20
import org.apache.logging.log4j.Logger;
21
import org.geotools.data.DataUtilities;
22
import org.geotools.feature.SchemaException;
23
import org.geotools.feature.simple.SimpleFeatureBuilder;
24
import org.geotools.geometry.jts.JTS;
25
import org.geotools.geometry.jts.JTSFactoryFinder;
26
import org.geotools.referencing.CRS;
27
import org.geotools.referencing.GeodeticCalculator;
28
import org.geotools.referencing.crs.DefaultGeographicCRS;
29
import org.locationtech.jts.geom.Coordinate;
30
import org.locationtech.jts.geom.CoordinateSequence;
31
import org.locationtech.jts.geom.Geometry;
32
import org.locationtech.jts.geom.GeometryFactory;
33
import org.locationtech.jts.geom.Point;
34
import org.locationtech.jts.geom.Polygon;
35
import org.locationtech.jts.util.GeometricShapeFactory;
36
import org.opengis.feature.GeometryAttribute;
37
import org.opengis.feature.simple.SimpleFeature;
38
import org.opengis.feature.simple.SimpleFeatureType;
39
import org.opengis.geometry.MismatchedDimensionException;
40
import org.opengis.referencing.FactoryException;
41
import org.opengis.referencing.NoSuchAuthorityCodeException;
42
import org.opengis.referencing.crs.CoordinateReferenceSystem;
43
import org.opengis.referencing.crs.ProjectedCRS;
44
import org.opengis.referencing.operation.MathTransform;
45
import org.opengis.referencing.operation.TransformException;
46

    
47
import si.uom.SI;
48
import tec.uom.se.quantity.Quantities;
49

    
50
/**
51
 * See :
52
 *
53
 * https://gis.stackexchange.com/questions/311272/create-dynamic-circle-polygon-from-specific-lat-long-using-geotools
54
 * https://gis.stackexchange.com/questions/283183/using-geometricshapefactory-to-create-circle-with-radius-in-miles
55
 *
56
 * https://stackoverflow.com/questions/36481651/how-do-i-create-a-circle-with-latitude-longitude-and-radius-with-geotools#36528805
57
 *
58
 *
59
 * Another great library for creating shapes is https://github.com/locationtech/spatial4j
60
 *
61
 * @author Andreas Kohlbecker
62
 * @since Apr 21, 2020
63
 */
64
public class GeometryBuilder {
65

    
66
	private final static Logger logger = LogManager.getLogger(GeometryBuilder.class);
67

    
68
	public enum CircleMethod {
69
		circle, simpleCircleSmall, simpleCircle, reprojectedCircle;
70
	}
71

    
72
	/**
73
	 * Fails with javax.measure.IncommensurableException: m is not compatible with
74
	 * deg
75
	 *
76
	 * see
77
	 * https://gis.stackexchange.com/questions/283183/using-geometricshapefactory-to-create-circle-with-radius-in-miles
78
	 *
79
	 * @param radius
80
	 * @param latitude
81
	 * @param longitude
82
	 * @return
83
	 * @throws NoSuchAuthorityCodeException
84
	 * @throws FactoryException
85
	 */
86
	public Polygon circle(Double radius, Double latitude, Double longitude)
87
			throws NoSuchAuthorityCodeException, FactoryException {
88

    
89
		Quantity<Length> radiusMeter = Quantities.getQuantity(radius, SI.METRE);
90
		Unit<Length> origUnit = (Unit<Length>) DefaultGeographicCRS.WGS84.getCoordinateSystem().getAxis(0).getUnit(); //// (Unit<Length>)
91
																														//// ////
92
																														//// origCRS.getCoordinateSystem().getAxis(0).getUnit();
93

    
94
		GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
95
		shapeFactory.setNumPoints(32);
96
		shapeFactory.setCentre(new Coordinate(latitude, longitude));
97
		shapeFactory.setSize(radiusMeter.to(origUnit).getValue().doubleValue() * 2);
98

    
99
		Polygon circlePolygon = shapeFactory.createCircle();
100
		circlePolygon.setSRID(4326);
101

    
102
		return circlePolygon;
103
	}
104

    
105
	/**
106
	 * Only suitable for small radius (> 1000 m) as the circles are heavily
107
	 * distorted otherwise.
108
	 *
109
	 * @param radius
110
	 * @param latitude
111
	 * @param longitude
112
	 * @return
113
	 */
114
	public Geometry simpleCircleSmall(Double radius, Double latitude, Double longitude) {
115

    
116
		double diameterInMeters = radius * 2.0d;
117
		GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
118
		shapeFactory.setNumPoints(64); // adjustable
119
		shapeFactory.setCentre(new Coordinate(longitude, latitude));
120
		// Length in meters of 1° of latitude = always 111.32 km
121
		shapeFactory.setWidth(diameterInMeters / 111320d);
122
		// Length in meters of 1° of longitude = 40075 km * cos( latitude ) / 360
123
		shapeFactory.setHeight(diameterInMeters / (40075000 * Math.cos(Math.toRadians(latitude)) / 360));
124

    
125
		Polygon circle = shapeFactory.createEllipse();
126
		return circle;
127
	}
128

    
129
	/**
130
	 * Creates perfect circles which are looking good but might be projected
131
	 * incorrectly for the resulting map
132
	 *
133
	 * @param distance
134
	 * @param latitude
135
	 * @param longitude
136
	 * @return
137
	 */
138
	public Geometry simpleCircle(Quantity<Length> distance, Double latitude, Double longitude) {
139

    
140
		GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
141
		CoordinateSequence coordinateSequence = geometryFactory.getCoordinateSequenceFactory()
142
				.create(new Coordinate[] { new Coordinate(longitude, latitude) });
143
		Point point = new Point(coordinateSequence, geometryFactory);
144

    
145
		GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
146
		calc.setStartingGeographicPoint(longitude, latitude);
147
		UnitConverter converter = distance.getUnit().getConverterTo(SI.METRE);
148
		double d = converter.convert(distance.getValue()).doubleValue();
149
		calc.setDirection(0.0, d);
150
		Point2D p2 = calc.getDestinationGeographicPoint();
151
		calc.setDirection(90.0, d);
152
		Point2D p3 = calc.getDestinationGeographicPoint();
153

    
154
		double dy = p2.getY() - latitude;
155
		double dx = p3.getX() - longitude;
156
		double dist = (dy + dx) / 2.0;
157
		Polygon p1 = (Polygon) point.buffer(dist);
158
		return p1;
159
	}
160

    
161
	/**
162
	 * This Method should produces the best circles.
163
	 *
164
	 * The code is based on an example published by Ian Turton on stackoverflow:
165
	 *
166
	 * https://stackoverflow.com/questions/36481651/how-do-i-create-a-circle-with-latitude-longitude-and-radius-with-geotools#36528805
167
	 *
168
	 * see https://gist.github.com/ianturton/973563fe5004985ba35a6e2247f7d823 and
169
	 * https://gitlab.com/snippets/17558
170
	 *
171
	 * @param feature
172
	 * @param distance
173
	 * @return
174
	 */
175
	public Geometry bufferFeature(SimpleFeature feature, Double radius) {
176

    
177
		// replacment for Measure<Double, Length> distance:
178
		Quantity<Length> radiusMeter = Quantities.getQuantity(radius, SI.METRE);
179

    
180
		// extract the geometry
181
		GeometryAttribute gProp = feature.getDefaultGeometryProperty();
182
		CoordinateReferenceSystem origCRS = gProp.getDescriptor().getCoordinateReferenceSystem();
183

    
184
		Geometry geom = (Geometry) feature.getDefaultGeometry();
185
		Geometry pGeom = geom;
186
		MathTransform toTransform, fromTransform = null;
187

    
188
		// reproject the geometry to a local projection
189
		if (!(origCRS instanceof ProjectedCRS)) {
190

    
191
			double x = geom.getCoordinate().x;
192
			double y = geom.getCoordinate().y;
193

    
194
			String code = "AUTO:42001," + x + "," + y;
195
			// System.out.println(code);
196
			CoordinateReferenceSystem auto;
197
			try {
198
				auto = CRS.decode(code);
199
				toTransform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);
200
				fromTransform = CRS.findMathTransform(auto, DefaultGeographicCRS.WGS84);
201
				pGeom = JTS.transform(geom, toTransform);
202

    
203
			} catch (MismatchedDimensionException | TransformException | FactoryException e) {
204
				logger.error(e);
205
			}
206

    
207
		}
208

    
209
		// create a buffer around the geometry, assumes the geometry is in the same
210
		// units as the distance variable.
211
		Geometry out = pGeom.buffer(radiusMeter.getValue().doubleValue());
212
		Geometry retGeom = out;
213
		// reproject the geometry to the original projection
214
		if (!(origCRS instanceof ProjectedCRS)) {
215
			try {
216
				retGeom = JTS.transform(out, fromTransform);
217

    
218
			} catch (MismatchedDimensionException | TransformException e) {
219
				logger.error(e);
220
			}
221
		}
222

    
223
		return retGeom;
224
	}
225

    
226
	public static SimpleFeature createSimplePointFeature(Double longitude, Double latitude) throws SchemaException {
227
		SimpleFeatureType schema = DataUtilities.createType("", "Location", "locations:Point:srid=4326,id:Integer" // a
228
																													// number
229
																													// attribute
230
		);
231
		SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(schema);
232

    
233
		GeometryFactory geometryFactory = new GeometryFactory();
234
		/* Longitude (= x coord) first ! */
235
		Point point = geometryFactory.createPoint(new Coordinate(longitude, latitude));
236
		featureBuilder.add(point);
237
		SimpleFeature feature = featureBuilder.buildFeature(null);
238
		return feature;
239
	}
240

    
241
}
(1-1/3)