Skip to content

Commit b16bec8

Browse files
Simone Giannecchinin-lagomarsini
Simone Giannecchini
authored and
n-lagomarsini
committed
[GEOT-4822] Modify RasterAsPointCollectionProcess for supporting new features.
1 parent 6c44f80 commit b16bec8

File tree

4 files changed

+1055
-116
lines changed

4 files changed

+1055
-116
lines changed

‎modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/GridConvergenceAngleCalc.java

Lines changed: 248 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,248 @@
1+
/*
2+
* GeoTools - The Open Source Java GIS Toolkit
3+
* http://geotools.org
4+
*
5+
* (C) 2014, Open Source Geospatial Foundation (OSGeo)
6+
* (C) 2014 TOPP - www.openplans.org.
7+
*
8+
* This library is free software; you can redistribute it and/or
9+
* modify it under the terms of the GNU Lesser General Public
10+
* License as published by the Free Software Foundation;
11+
* version 2.1 of the License.
12+
*
13+
* This library is distributed in the hope that it will be useful,
14+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16+
* Lesser General Public License for more details.
17+
*/
18+
package org.geotools.process.raster;
19+
20+
import org.geotools.geometry.DirectPosition2D;
21+
import org.geotools.referencing.GeodeticCalculator;
22+
import org.opengis.geometry.DirectPosition;
23+
import org.opengis.referencing.crs.CoordinateReferenceSystem;
24+
import org.opengis.referencing.cs.AxisDirection;
25+
import org.opengis.referencing.cs.CoordinateSystemAxis;
26+
27+
/**
28+
* Used to calculate an estimate for the Grid Convergence Angle at a point within a 2D Coordinate Reference System. The Grid Convergence Angle at any
29+
* point on a projected 2D map is the difference between "up" or "grid north" on the map at that point and True North. Since the map is projected, the
30+
* Grid Convergence Angle can change at each point on the map; True North can be in a different Cartesian direction on the flat map for every point.
31+
* One example impact of this is that vectors cannot be accurately drawn on the screen by simply rotating them a certain amount in "screen degrees."
32+
* This class, then, is used to estimate the Grid Convergence Angle at those points. Though the underlying meaning is the same, different mapping
33+
* conventions define the angle's direction (+/-) and beginning point / ending point differently. Therefore, for the purposes of this code, the Grid
34+
* Convergence Angle is defined as the angle (C) FROM true north TO grid north with 0 deg up and angles increasing positively clockwise. So, for the
35+
* example below, C would be approximately -33 deg, since the angle FROM true north TO grid north is Counter-Clockwise.
36+
*
37+
* <pre>
38+
*
39+
* Up
40+
* Grid North
41+
* 0 deg True North
42+
* | .
43+
* | C .
44+
* | .
45+
* | .
46+
* 270 deg-----------O--------------90 deg
47+
* |
48+
* |
49+
* |
50+
* |
51+
* |
52+
* 180 deg
53+
*
54+
* </pre>
55+
*
56+
* This class uses the GeoTools GeodeticCalculator for performing underlying calculations. <br>
57+
* <br>
58+
* Some literature (don't have a link) says that the Grid Covergence Angle is really the angle between a line extending toward grid north and one
59+
* extending toward true north THAT HAVE BEEN PROJECTED into the map projection, but suggests the difference between this angle and the way the angle
60+
* is being estimated here is small. May want some verification of that. <br>
61+
* <br>
62+
*
63+
* @author Mike Grogan, WeatherFlow, Inc., Synoptic <br>
64+
* <br>
65+
* @see http://www.threelittlemaids.co.uk/magdec/explain.html
66+
* @see http://www.bluemarblegeo.com/knowledgebase/calculator/Scale_Factor_and_Convergence.htm
67+
*
68+
*/
69+
final class GridConvergenceAngleCalc {
70+
71+
/** Input Coverage CRS */
72+
private final CoordinateReferenceSystem crs;
73+
74+
/** Operator used for calculating the GridConvergence angle */
75+
private GeodeticCalculator geoCalc;
76+
77+
/** Index associated to the "up" axis */
78+
private final int upAxisDimension;
79+
80+
/**
81+
* Constructs a new GridConvergenceAngleCalc for a given CoordinateReferenceSystem
82+
*
83+
* @param crs CoordinateReferenceSystem for which to construct a new GridConvergenceAngleCalc.
84+
* @throws Exception
85+
*/
86+
public GridConvergenceAngleCalc(CoordinateReferenceSystem crs) throws Exception {
87+
this.crs = crs;
88+
this.geoCalc = new GeodeticCalculator(this.crs);
89+
this.upAxisDimension = determineUpAxisDimension();
90+
//
91+
// If we could not find the "up" axis ... meaning up on the map/screen
92+
// not in the vertical ... then throw an exception
93+
//
94+
95+
if (upAxisDimension < 0) {
96+
throw new Exception("Up Axis can not be determined.");
97+
}
98+
}
99+
100+
/**
101+
* Estimates the grid convergence angle at a position within a Coordinate Reference System. The angle returned is as described in the
102+
* documentation for the class. The Coordinate Reference System of the supplied position must be the same as was used when constructing the
103+
* calculator, because using anything else would not make sense as convergence angle depends on projection.
104+
*
105+
* @param position DirectPosition2D at which we want to estimate the grid convergence angle
106+
* @return double containing grid convergence angle, as described in documentation for the class.
107+
* @throws Exception
108+
*/
109+
public double getConvergenceAngle(DirectPosition2D position) throws Exception {
110+
111+
//
112+
// Check to make sure the coordinate reference system for the
113+
// argument is the same as the calculator.
114+
//
115+
116+
CoordinateReferenceSystem positionCRS = position.getCoordinateReferenceSystem();
117+
118+
if (!positionCRS.equals(crs)) {
119+
throw new Exception("Position CRS does not match Calculator CRS");
120+
}
121+
122+
//
123+
// We will use the Geotools Geodetic calculator to estimate the
124+
// convergence angle. We estimate this by taking the supplied point,
125+
// moving "upward" along the proper upward map axis by 1 unit, and
126+
// then having the Geodetic calculator tell us the azimuth from the
127+
// starting point to the ending point. Since the azimuth is relative
128+
// to true north ... and we are "walking" along a grid north
129+
// parallel, the azimuth then essentially tells us the angle from
130+
// true north to grid north, or the local grid convergence angle.
131+
//
132+
133+
//
134+
// Get the "up" axis
135+
//
136+
137+
CoordinateSystemAxis upAxis = crs.getCoordinateSystem().getAxis(upAxisDimension);
138+
139+
//
140+
// Need to make sure we're not going to go out of bounds along the
141+
// axis by going up a little bit.
142+
//
143+
// Determine the maximum value along that axis
144+
//
145+
146+
double upAxisMax = upAxis.getMaximumValue();
147+
148+
//
149+
// Get the starting value along the up axis
150+
//
151+
152+
double startValueUp = position.getOrdinate(upAxisDimension);
153+
154+
//
155+
// If adding 1 to the up axis is going to push us out of bounds, then
156+
// first subtract 1 from the starting position ... the estimate should
157+
// still be close if units are close.
158+
//
159+
160+
if ((startValueUp + 1) > upAxisMax) {
161+
position.setOrdinate(upAxisDimension, position.getOrdinate(upAxisDimension) - 1);
162+
}
163+
164+
//
165+
// Set the starting position for the geodetic calculator to position.
166+
//
167+
168+
geoCalc.setStartingPosition(position);
169+
170+
//
171+
// Set the ending position to be the same as the starting position,
172+
// except move "up" 1 unit along the "up" axis.
173+
//
174+
175+
DirectPosition2D endingPosition = new DirectPosition2D((DirectPosition) position);
176+
endingPosition.setOrdinate(upAxisDimension, position.getOrdinate(upAxisDimension) + 1);
177+
geoCalc.setDestinationPosition(endingPosition);
178+
179+
//
180+
// Now just ask for the azimuth, which is our convergence angle
181+
// estimate.
182+
//
183+
return geoCalc.getAzimuth();
184+
}
185+
186+
/**
187+
* Determines which axis in the calculator's Coordinate Reference System is "up"
188+
*
189+
* @return int with up axis dimension, or -1 if up axis cannot be found.
190+
*/
191+
private int determineUpAxisDimension() {
192+
//
193+
// Grab the number of dimensions. We only can deal with a 2D
194+
// projection here. Set to -1 if not a 2D system ... and let
195+
// other code throw errors.
196+
//
197+
198+
int numDimensions = crs.getCoordinateSystem().getDimension();
199+
200+
if (numDimensions > 2) {
201+
return -1;
202+
}
203+
204+
//
205+
// Loop through all of the axes until you find the one that is
206+
// probably the upward axis.
207+
//
208+
209+
for (int i = 0; i < numDimensions; i++) {
210+
CoordinateSystemAxis axis = crs.getCoordinateSystem().getAxis(i);
211+
AxisDirection axisDirection = axis.getDirection();
212+
213+
if (axisDirection.equals(AxisDirection.DISPLAY_UP)
214+
|| axisDirection.equals(AxisDirection.EAST_NORTH_EAST)
215+
|| axisDirection.equals(AxisDirection.NORTH)
216+
|| axisDirection.equals(AxisDirection.NORTH_EAST)
217+
|| axisDirection.equals(AxisDirection.NORTH_NORTH_EAST)
218+
|| axisDirection.equals(AxisDirection.NORTH_NORTH_WEST)
219+
|| axisDirection.equals(AxisDirection.NORTH_WEST)
220+
|| axisDirection.equals(AxisDirection.ROW_POSITIVE)
221+
|| axisDirection.equals(AxisDirection.UP)
222+
|| axisDirection.equals(AxisDirection.WEST_NORTH_WEST)) {
223+
return i;
224+
}
225+
}
226+
227+
//
228+
// If not found yet, find one with name Northing or y and assume
229+
// it is up.
230+
//
231+
232+
for (int i = 0; i < numDimensions; i++) {
233+
CoordinateSystemAxis axis = crs.getCoordinateSystem().getAxis(i);
234+
String axisName = axis.getName().toString().toUpperCase();
235+
if (axisName.equals("Y") || axisName.equals("NORTHING")
236+
|| axisName.contains("NORTHING")) {
237+
return i;
238+
}
239+
}
240+
241+
//
242+
// If the up axis still hasn't been found, then signify we can't
243+
// find it with a -1.
244+
//
245+
246+
return -1;
247+
}
248+
}

0 commit comments

Comments
 (0)