Schlagwort-Archive: Projection

GIS in RapidMiner (3) – Distanz, Fläche

(English version)

Nach der Einführung und dem Datenimport geht es nun an echte geographische Berechnungen.

Eine der wichtigsten Informationen ist die Distanz von Objekten voneinander oder einem Referenzpunkt. Auch Data-Mining-Verfahren wie k Nearest Neighbors berechnen Distanzen.

Zwei verschiedene Methoden stehen uns zur Verfügung, um Distanzen zwischen zwei Punkten auf der Erdoberfläche zu berechnen: Entweder können wir die Punkte als zweidimensionale Geometrie mit X- und Y-Koordinaten auffassen (die Berechnung ist dann ganz einfach), oder die Distanz auf der Oberfläche des Ellipsoids der Erde ausrechnen. Die zweite Vorgehensweise ist mathematisch natürlich deutlich aufwändiger, liefert aber bei größeren Distanzen (z. B. Orte auf verschiedenen Kontinenten) genauere Daten. Deswegen wird in vielen Anwendungen, die nur Objekte in einem eingeschränkten Gebiet enthalten, auf die erste Methode zurückgegriffen.

Transformation in eine andere Projektion (Koordinatensystem)

Wie in der Einführung ausgeführt müssen wir die Geometrien manchmal in andere Projektionen transformieren, um mit sinnvollen Einheiten wie Meter rechnen zu können. Die Methoden dafür sind in GeoScript enthalten und ihre Anwendung ist recht einfach:


import geoscript.proj.*;

fromProj = new Projection("epsg:4326");
toProj = new Projection("epsg:3035");

projectedGeom = Projection.transform(geom, fromProj, toProj);

Ein fertig anwendbarer RapidMiner-Prozess befindet sich hier. Er braucht drei Parameter, die als Makros im Prozesskontext definiert sind und beim Aufruf angegeben werden können:

GEOM_ATT: Name des Attributs, das die zu transformierende Geometrie (im WKT-Format) enthält

FROM_PROJECTION, TO_PROJECTION: Die EPSG-Nummern der Quell- und Zielprojektion.

Damit lassen sich die Koordinaten leicht von einem allgemeinen Koordinatensystem wie WGS84 (Längen- und Breitengrad, EPSG:4326) in  ein gebräuchlicheres wie z. B. ETRS89/Austria Lambert (EPSG:3416) konvertieren. Das werden wir im nächsten Beispielprozess anwenden, um Distanzen zwischen Objekten in Wien, aber auch Fläche und Umfang von Bezirken zu berechnen.

Berechnung von Distanz, Fläche und Umfang

GeoScript enthält dafür einfach anzuwendende Funktionen:

flaeche = geom.getArea();
umfang = geom.getLength();
//Für die Distanz brauchen wir eine zweite Geometrie
distanz = geom1.distance(geom2);

Sobald man das Geometrie-Objekt in einer richtigen Projektion hat, sind die Berechnungen ganz simpel.

getArea() liefert die Fläche eines Polygons oder einer Polygongruppe (Multipolygon); getLength() die Länge einer Linie oder den Umfang eines Polygons, jeweils in den Einheiten der aktuellen Projektion.

Für die Distanz brauchen wir ein zweites Geometrie-Objekt, das nicht unbedingt ein Punkt sein muß – es ist auch möglich, die Distanz zwischen Linien und Flächen zu berechnen.

Der Beispielprozess holt zwei Datensätze vom Open-Data-Portal der Stadt Wien: Museen (Punkte) und Bezirksgrenzen (Polygone). Beide werden aus Längen- und Breitengrad-Koordinaten in eine in Österreich gebräuchliche, meter-basierte Projektion transformiert. Mit getArea() und getLength() werden Fläche und Umfang der Bezirke berechnet. Da der Original-Datensatz diese Information bereits enthält, können wir leicht prüfen, ob die Berechnungen korrekt sind. (Kleine Unterschiede resultieren wohl aus der Rundung der Koordinaten für den CSV-Export.)

Dann wird noch der erste Bezirk selektiert und mit dem Museums-Datensatz zusammengeführt. In diesem kombinierten Datensatz haben wir nun zwei Geometrien, wir können also die Entfernung des Museums von der Innenstadt berechnen. Punkte, die im Polygon des ersten Bezirks liegen, haben die Distanz 0.

Berechnung von Distanzen auf der Erdoberfläche

Die zweite, genauere, aber langsamere Berechnungsmethode kann auch recht einfach angewendet werden. Hierfür importieren wir aus der GeoTools-Library, auf die GeoScript aufbaut, die Klasse GeodeticCalculator. (Die Bibliotheken, die wir in der Einführung für GeoScript übernommen haben, reichen dafür aus, wir müssen also nichts Neues installieren.)

Für diese Methode brauchen wir die Längen- und Breitengrade von zwei Punkten, also WGS84-Koordinaten. (Es gäbe auch die Möglichkeit, transformierte Koordinaten zu verwenden, dafür müßten wir dem GeodeticCalculator auch das Koordinatensystem übergeben.)

import org.geotools.referencing.GeodeticCalculator;

gcalc = new GeodeticCalculator();

gcalc.setStartingGeographicPoint(lon1, lat1);
gcalc.setDestinationGeographicPoint(lon2, lat2);

distance = gcalc.getOrthodromicDistance();

Hier ist es wichtig, die Reihenfolge der Koordinatenangaben zu beachten. In anderen Bereichen sind wir gewohnt, X- und Y-Koordinaten in dieser Reihenfolge anzugeben, GIS-Systeme arbeiten jedoch manchmal mit der Reihenfolge Y, X.

Der Beispielprozess enthält die Koordinatenpaare einiger Hauptstädte auf verschiedenen Kontinenten. Mit einem Cartesian Product werden alle Städte mit allen anderen verknüpft und jeweils die Distanzen in km berechnet. (Die berechneten Distanzen habe ich mit PostGIS verifiziert; die Ergebnisse sind sehr genau.) Für diese Distanzen würde eine Berechnung mit der Geometrie-Methode schon sehr große Ungenauigkeiten liefern, hier empfiehlt es sich also sehr, die GeodeticCalculator-Methode zu nutzen.

GIS in RapidMiner (3) – Distance and Area

After the introduction and data import we can start to perform actual geographic calculations.

One of the most important measures is the distance of objects from each other or from a reference point. Distances are also calculated by data mining operations like k Nearest Neighbors.

Two different methods of calculating distances between points on the surface of Earth exist: Either we can pretend that the points are in a two-dimensional geometry with X and Y coordinates (which makes the distance calculation very easy), or we use the actual Earth ellipsoid for the calculation. The second method is very computing-intensive, but it returns more precise results when used on larger distances (e. g. places on different continents). For many operations acting on a limited area, the first method is used.

Transformation to another projection (coordinate system)

As described in the introduction, we often need to transform geometries to different projections so we can use units like meters. Coordinate system transformation is also available in GeoScript, and using it is quite easy.

import geoscript.proj.*;

//Source and destination projections
fromProj = new Projection("epsg:4326");
toProj = new Projection("epsg:3035");

projectedGeom = Projection.transform(geom, fromProj, toProj);

Here is a readily usable RapidMiner process. It takes three parameters that can be specified as macros in the process context. You can overwrite them when calling the process in the Execute Process operator.

GEOM_ATT: Name of the attribute that contains the geometry to be transformed (in WKT format)

FROM_PROJECTION, TO_PROJECTION: The EPSG numbers of the source and target projections

With this process, you can easily transform geometries from a common coordinate system like WGS84 (Latitude and Longitude, EPSG:4326) to a special one, for example ETRS89/Austria Lambert (EPSG:3416). We will use this in the example process for calculating distances between objects in Vienna, Austria and also determine the area and circumference of Vienna’s 23 districts.

Calculating distance, area and circumference

GeoScript contains easy to use functions:

area = geom.getArea();
circumference = geom.getLength();
//We need a second geometry for calculating distance
distance = geom1.distance(geom2);

After having transformed the geometry object into a matching projection, the calculations are really simple.

getArea() returns the area of a polygon or a group of polygons (Multipolygon); getLength() gives the length of a line or the circumference of a polygon in the units of the used projection.

For calculating the distance, a second geometry object is required. It doesn’t need to be a point: it’s also possible to calculate the distance between lines and areas.

The example process fetches two data sets from the Vienna Open Data Portal: Museums (points) and district borders (polygons). The process transforms then from latitude and longitude coordinates into a projection used in Austria. One script calculates the area and circumference of the districts using getArea() and getLength(). The original data set already contains these measures so we can easily check that they’re correct. (Small differences are the consequence of rounding the coordinates for CSV export.)

After that, the first district (Inner City) is selected and joined with the museum data set. The combined data set contains two geometries in a common projection, so we can calculate the distance between the museum and the Inner City. Points in the first district have a distance of 0.

Calculating distances on the surface of Earth

The second calculation method (more precise but slower) can also be used quite easily. For this, we import the class GeodeticCalculator from the GeoTools library, a base component of GeoScript. (The GeoScript libraries installed in the introduction are enough for this, we don’t need to set up more stuff.)

For using this method, we need latitude and longitudes of two points, in other words WGS84 coordinates. (It would be possible to use transformed coordinates by specifying a coordinate system in the GeodeticCalculator.)

import org.geotools.referencing.GeodeticCalculator;

gcalc = new GeodeticCalculator();

gcalc.setStartingGeographicPoint(lon1, lat1);
gcalc.setDestinationGeographicPoint(lon2, lat2);

distance = gcalc.getOrthodromicDistance();

Be careful when specifying the coordinates. We usually write X and Y coordinates in this order but GIS tools often use the order Y, X.

The example process contains coordinate pairs of a few capital cities lying on different continents. We build a Cartesian Product of all cities and calculate their distances in kilometers. (The distances were verified with PostGIS, they are very precise.) On these distances, using the geometry method would result in huge inaccuracies, so it’s really better to use the GeodeticCalculator method there.