ashatora/openlayers
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Adding methods for getting geodesic measures from geometries. Assuming geometries can be transformed into Geographic/WGS84, getGeodesicLength and getGeodesicArea should return reasonable 'on the ground' metrics. Use getLength and getArea for the planar metrics. r=crschmidt (closes #1819)

Browse Source 
git-svn-id: http://svn.openlayers.org/trunk/openlayers@9248 dc9f47b5-9b13-0410-9fdd-eb0c1a62fdaf
This commit is contained in:
Tim Schaub
2009-04-08 23:12:24 +00:00
parent 5f335f4207
commit 08077f0f42
8 changed files with 250 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

63
examples/measure.html
View File

@@ -64,12 +64,13 @@
measureControls = { measureControls = {
line: new OpenLayers.Control.Measure( line: new OpenLayers.Control.Measure(
OpenLayers.Handler.Path, { OpenLayers.Handler.Path, {
persist: true, persist: true,
handlerOptions: { handlerOptions: {
layerOptions: {styleMap: styleMap}, layerOptions: {styleMap: styleMap}
}
} }
}), ),
polygon: new OpenLayers.Control.Measure( polygon: new OpenLayers.Control.Measure(
OpenLayers.Handler.Polygon, { OpenLayers.Handler.Polygon, {
persist: true, persist: true,
@@ -94,31 +95,6 @@
document.getElementById('noneToggle').checked = true; document.getElementById('noneToggle').checked = true;
} }
function calcVincenty(geometry) {
/**
* Note: this function assumes geographic coordinates and
* will fail otherwise. OpenLayers.Util.distVincenty takes
* two objects representing points with geographic coordinates
* and returns the geodesic distance between them (shortest
* distance between the two points on an ellipsoid) in *kilometers*.
*
* It is important to realize that the segments drawn on the map
* are *not* geodesics (or "great circle" segments). This means
* that in general, the measure returned by this function
* will not represent the length of segments drawn on the map.
*/
var dist = 0;
for (var i = 1; i < geometry.components.length; i++) {
var first = geometry.components[i-1];
var second = geometry.components[i];
dist += OpenLayers.Util.distVincenty(
{lon: first.x, lat: first.y},
{lon: second.x, lat: second.y}
);
}
return dist;
}
function handleMeasurements(event) { function handleMeasurements(event) {
var geometry = event.geometry; var geometry = event.geometry;
@@ -129,10 +105,6 @@
var out = ""; var out = "";
if(order == 1) { if(order == 1) {
out += "measure: " + measure.toFixed(3) + " " + units; out += "measure: " + measure.toFixed(3) + " " + units;
if (map.getProjection() == "EPSG:4326") {
out += "<br /> Great Circle Distance: " +
calcVincenty(geometry).toFixed(3) + " km *";
}
} else { } else {
out += "measure: " + measure.toFixed(3) + " " + units + "<sup>2</" + "sup>"; out += "measure: " + measure.toFixed(3) + " " + units + "<sup>2</" + "sup>";
} }
@@ -149,6 +121,13 @@
} }
} }
} }
function toggleGeodesic(element) {
for(key in measureControls) {
var control = measureControls[key];
control.geodesic = element.checked;
}
}
</script> </script>
</head> </head>
<body onload="init()"> <body onload="init()">
@@ -174,13 +153,17 @@
<input type="radio" name="type" value="polygon" id="polygonToggle" onclick="toggleControl(this);" /> <input type="radio" name="type" value="polygon" id="polygonToggle" onclick="toggleControl(this);" />
<label for="polygonToggle">measure area</label> <label for="polygonToggle">measure area</label>
</li> </li>
<li>
<input type="checkbox" name="geodesic" id="geodesicToggle" onclick="toggleGeodesic(this);" />
<label for="polygonToggle">use geodesic measures</label>
</li>
</ul> </ul>
<p>* Note that the geometries drawn are planar geometries and the <p>Note that the geometries drawn are planar geometries and the
metrics returned by the measure control are planar measures. The metrics returned by the measure control are planar measures by
"great circle" distance does not necessarily represent the length default. If your map is in a geographic projection or you have the
of the segments drawn on the map. Instead, it is a geodesic metric that appropriate projection definitions to transform your geometries into
represents the cumulative shortest path between all vertices in the geographic coordinates, you can set the "geodesic" property of the control
geometry were they projected onto a sphere.</p> to true to calculate geodesic measures instead of planar measures.</p>
</div> </div>
</body> </body>
</html> </html>

30
lib/OpenLayers/Control/Measure.js
View File

@@ -56,6 +56,14 @@ OpenLayers.Control.Measure = OpenLayers.Class(OpenLayers.Control, {
*/ */
displaySystem: 'metric', displaySystem: 'metric',
/**
* Property: geodesic
* {Boolean} Calculate geodesic metrics instead of planar metrics. This
* requires that geometries can be transformed into Geographic/WGS84
* (if that is not already the map projection). Default is false.
*/
geodesic: false,
/** /**
* Property: displaySystemUnits * Property: displaySystemUnits
* {Object} Units for various measurement systems. Values are arrays * {Object} Units for various measurement systems. Values are arrays
@@ -213,10 +221,17 @@ OpenLayers.Control.Measure = OpenLayers.Class(OpenLayers.Control, {
* {Float} The geometry area in the given units. * {Float} The geometry area in the given units.
*/ */
getArea: function(geometry, units) { getArea: function(geometry, units) {
var area = geometry.getArea(); var area, geomUnits;
if(this.geodesic) {
area = geometry.getGeodesicArea(this.map.getProjectionObject());
geomUnits = "m";
} else {
area = geometry.getArea();
geomUnits = this.map.getUnits();
}
var inPerDisplayUnit = OpenLayers.INCHES_PER_UNIT[units]; var inPerDisplayUnit = OpenLayers.INCHES_PER_UNIT[units];
if(inPerDisplayUnit) { if(inPerDisplayUnit) {
var inPerMapUnit = OpenLayers.INCHES_PER_UNIT[this.map.getUnits()]; var inPerMapUnit = OpenLayers.INCHES_PER_UNIT[geomUnits];
area *= Math.pow((inPerMapUnit / inPerDisplayUnit), 2); area *= Math.pow((inPerMapUnit / inPerDisplayUnit), 2);
} }
return area; return area;
@@ -257,10 +272,17 @@ OpenLayers.Control.Measure = OpenLayers.Class(OpenLayers.Control, {
* {Float} The geometry length in the given units. * {Float} The geometry length in the given units.
*/ */
getLength: function(geometry, units) { getLength: function(geometry, units) {
var length = geometry.getLength(); var length, geomUnits;
if(this.geodesic) {
length = geometry.getGeodesicLength(this.map.getProjectionObject());
geomUnits = "m";
} else {
length = geometry.getLength();
geomUnits = this.map.getUnits();
}
var inPerDisplayUnit = OpenLayers.INCHES_PER_UNIT[units]; var inPerDisplayUnit = OpenLayers.INCHES_PER_UNIT[units];
if(inPerDisplayUnit) { if(inPerDisplayUnit) {
var inPerMapUnit = OpenLayers.INCHES_PER_UNIT[this.map.getUnits()]; var inPerMapUnit = OpenLayers.INCHES_PER_UNIT[geomUnits];
length *= (inPerMapUnit / inPerDisplayUnit); length *= (inPerMapUnit / inPerDisplayUnit);
} }
return length; return length;

46
lib/OpenLayers/Geometry/Collection.js
View File

@@ -231,6 +231,52 @@ OpenLayers.Geometry.Collection = OpenLayers.Class(OpenLayers.Geometry, {
return area; return area;
}, },
/**
* APIMethod: getGeodesicArea
* Calculate the approximate area of the polygon were it projected onto
* the earth.
*
* Parameters:
* projection - {<OpenLayers.Projection>} The spatial reference system
* for the geometry coordinates. If not provided, Geographic/WGS84 is
* assumed.
*
* Reference:
* Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
* Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
* Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
*
* Returns:
* {float} The approximate geodesic area of the geometry in square meters.
*/
getGeodesicArea: function(projection) {
var area = 0.0;
for(var i=0, len=this.components.length; i<len; i++) {
area += this.components[i].getGeodesicArea(projection);
}
return area;
},
/**
* APIMethod: getGeodesicLength
* Calculate the approximate length of the geometry were it projected onto
* the earth.
*
* projection - {<OpenLayers.Projection>} The spatial reference system
* for the geometry coordinates. If not provided, Geographic/WGS84 is
* assumed.
*
* Returns:
* {Float} The appoximate geodesic length of the geometry in meters.
*/
getGeodesicLength: function(projection) {
var length = 0.0;
for(var i=0, len=this.components.length; i<len; i++) {
length += this.components[i].getGeodesicLength(projection);
}
return length;
},
/** /**
* APIMethod: move * APIMethod: move
* Moves a geometry by the given displacement along positive x and y axes. * Moves a geometry by the given displacement along positive x and y axes.

36
lib/OpenLayers/Geometry/Curve.js
View File

@@ -52,5 +52,41 @@ OpenLayers.Geometry.Curve = OpenLayers.Class(OpenLayers.Geometry.MultiPoint, {
return length; return length;
}, },
/**
* APIMethod: getGeodesicLength
* Calculate the approximate length of the geometry were it projected onto
* the earth.
*
* projection - {<OpenLayers.Projection>} The spatial reference system
* for the geometry coordinates. If not provided, Geographic/WGS84 is
* assumed.
*
* Returns:
* {Float} The appoximate geodesic length of the geometry in meters.
*/
getGeodesicLength: function(projection) {
var geom = this; // so we can work with a clone if needed
if(projection) {
var gg = new OpenLayers.Projection("EPSG:4326");
if(!gg.equals(projection)) {
geom = this.clone().transform(projection, gg);
}
}
var length = 0.0;
if(geom.components && (geom.components.length > 1)) {
var p1, p2;
for(var i=1, len=geom.components.length; i<len; i++) {
p1 = geom.components[i-1];
p2 = geom.components[i];
// this returns km and requires lon/lat properties
length += OpenLayers.Util.distVincenty(
{lon: p1.x, lat: p1.y}, {lon: p2.x, lat: p2.y}
);
}
}
// convert to m
return length * 1000;
},
CLASS_NAME: "OpenLayers.Geometry.Curve" CLASS_NAME: "OpenLayers.Geometry.Curve"
}); });

44
lib/OpenLayers/Geometry/LinearRing.js
View File

@@ -205,6 +205,50 @@ OpenLayers.Geometry.LinearRing = OpenLayers.Class(
return area; return area;
}, },
/**
* APIMethod: getGeodesicArea
* Calculate the approximate area of the polygon were it projected onto
* the earth. Note that this area will be positive if ring is oriented
* clockwise, otherwise it will be negative.
*
* Parameters:
* projection - {<OpenLayers.Projection>} The spatial reference system
* for the geometry coordinates. If not provided, Geographic/WGS84 is
* assumed.
*
* Reference:
* Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
* Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
* Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
*
* Returns:
* {float} The approximate signed geodesic area of the polygon in square
* meters.
*/
getGeodesicArea: function(projection) {
var ring = this; // so we can work with a clone if needed
if(projection) {
var gg = new OpenLayers.Projection("EPSG:4326");
if(!gg.equals(projection)) {
ring = this.clone().transform(projection, gg);
}
}
var area = 0.0;
var len = ring.components && ring.components.length;
if(len > 2) {
var p1, p2;
for(var i=0; i<len-1; i++) {
p1 = ring.components[i];
p2 = ring.components[i+1];
area += OpenLayers.Util.rad(p2.x - p1.x) *
(2 + Math.sin(OpenLayers.Util.rad(p1.y)) +
Math.sin(OpenLayers.Util.rad(p2.y)));
}
area = area * 6378137.0 * 6378137.0 / 2.0;
}
return area;
},
/** /**
* Method: containsPoint * Method: containsPoint
* Test if a point is inside a linear ring. For the case where a point * Test if a point is inside a linear ring. For the case where a point

29
lib/OpenLayers/Geometry/Polygon.js
View File

@@ -60,6 +60,35 @@ OpenLayers.Geometry.Polygon = OpenLayers.Class(
return area; return area;
}, },
/**
* APIMethod: getGeodesicArea
* Calculate the approximate area of the polygon were it projected onto
* the earth.
*
* Parameters:
* projection - {<OpenLayers.Projection>} The spatial reference system
* for the geometry coordinates. If not provided, Geographic/WGS84 is
* assumed.
*
* Reference:
* Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
* Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
* Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
*
* Returns:
* {float} The approximate geodesic area of the polygon in square meters.
*/
getGeodesicArea: function(projection) {
var area = 0.0;
if(this.components && (this.components.length > 0)) {
area += Math.abs(this.components[0].getGeodesicArea(projection));
for(var i=1, len=this.components.length; i<len; i++) {
area -= Math.abs(this.components[i].getGeodesicArea(projection));
}
}
return area;
},
/** /**
* Method: containsPoint * Method: containsPoint
* Test if a point is inside a polygon. Points on a polygon edge are * Test if a point is inside a polygon. Points on a polygon edge are

23
tests/Geometry/LineString.html
View File

@@ -351,7 +351,30 @@
"clone() creates an OpenLayers.Geometry.LineString"); "clone() creates an OpenLayers.Geometry.LineString");
t.ok(geometry.equals(clone), "clone has equivalent coordinates"); t.ok(geometry.equals(clone), "clone has equivalent coordinates");
} }
function test_getGeodesicLength(t) {
// expected values from http://www.movable-type.co.uk/scripts/latlong-vincenty.html
var cases = [{
wkt: "LINESTRING(0 0, -10 45)",
exp: 5081689.690
}, {
wkt: "LINESTRING(-10 45, 0 0)",
exp: 5081689.690
}, {
wkt: "LINESTRING(0 0, -10 45, -20 50)",
exp: 5081689.690 + 935018.062
}];
t.plan(cases.length);
var geom, got;
for(var i=0; i<cases.length; ++i) {
geom = new OpenLayers.Geometry.fromWKT(cases[i].wkt);
got = geom.getGeodesicLength();
t.eq(Math.round(got), Math.round(cases[i].exp), "[case " + i + "] length calculated");
}
}
</script> </script>
</head> </head>

23
tests/Geometry/Polygon.html
View File

File diff suppressed because one or more lines are too long