Archiv für ‘Python’

September 24th, 2013

pyMCFsimplex – Efficient Solving of Minimum Cost Flow Problems in Python

pyMCFsimplex – a Python Wrapper for MCFSimplex

pyMCFimplex is a Python-Wrapper for the C++ MCFSimplex Solver Class from the Operations Research Group at the University of Pisa. MCFSimplex is a piece of software hat solves big sized Minimum Cost Flow Problems very fast through the (primal or dual) network simplex algorithm. Also important: MCFSimplex is open source software – it has been released under the LGPL license.

To get a feeling for the performance of MCFSimplex: it solves a MCFP of 50.000 nodes and 500.000 arcs in 5 seconds on a average linux box (2,2Ghz 6GB RAM, Ubuntu 64Bit). Another one: NetworkX (a pure python graph library) also has an implementation of the network simplex algorithm – but this one is magnitudes slower than e.g. MCFSimplex. I tried so solve the same MCFP (n = 50.000, m = 500.000) with NetworkX and terminated the process after 1 hour. I intend to publish more on performance tests in an upcoming article.
See also this paper for a MCFP benchmarks.

You may ask: So what about pyMCFsimplex? Well, you’ll get the performance of a well coded C++ MCFP solver with the simplicity of a Python API. Isn’t this nice?

weiterlesen »

2 Leute mögen diesen Artikel.

April 10th, 2013

Solving the Minimum Cost Flow Problem (6) – Google or-tools

Push-relabel algorithm: Google or-tools

At last we will solve the instance of a Minimum cost flow problem described in (1) with Google or-tools. Node 1 is the source node, nodes 2 and 3 are the transshipment nodes and node 4 is the sink node.

mcf_dimacs_1

So let’s see how the python API of google or-tools work:

# import graph from google or-tools
from graph import pywrapgraph

# Create graph
# Attention: add +1 for number of nodes and arcs, if your node-IDs begin with 1
# because google or-tools mcf solver internal counters are strictly zero-based!
num_nodes = 4 + 1
num_arcs = 5 + 1
# args: NumNodes * NumArcs
G = pywrapgraph.StarGraph(num_nodes, num_arcs)

# create min cost flow solver 
min_cost_flow = pywrapgraph.MinCostFlow(G)

# add node to graph with positive supply for each supply node 
min_cost_flow.SetNodeSupply(1, 4)
# add node to graph with negative demand for each demand node 
min_cost_flow.SetNodeSupply(4, -4)
# you can ignore transshipment nodes with zero supply when you are working with
# the mcfp-solver of google or-tools

# add arcs to the graph
arc = G.AddArc(1, 2)
min_cost_flow.SetArcCapacity(arc, 4)
min_cost_flow.SetArcUnitCost(arc, 2)
arc = G.AddArc(1, 3)
min_cost_flow.SetArcCapacity(arc, 2)
min_cost_flow.SetArcUnitCost(arc, 2)
arc = G.AddArc(2, 3)
min_cost_flow.SetArcCapacity(arc, 2)
min_cost_flow.SetArcUnitCost(arc, 1)
arc = G.AddArc(2, 4)
min_cost_flow.SetArcCapacity(arc, 3)
min_cost_flow.SetArcUnitCost(arc, 3)
arc = G.AddArc(3, 4)
min_cost_flow.SetArcCapacity(arc, 5)
min_cost_flow.SetArcUnitCost(arc, 1)

# solve the min cost flow problem
# flowDict contains the optimized flow
# flowCost contains the total minimized optimum
min_cost_flow.Solve()
flowCost = min_cost_flow.GetOptimalCost()

print "Optimum: %s" %flowCost

As one can see the python API of google or-tools is also pretty handy although everything is coded in C++. The python API was created using SWIG. The biggest caveat one can run into is the fact that the MCFP solver in google or-tools is strictly zero based. So adding nodes starting from 1 ..n result in strange errors if you don’t size the initial graph properly: just add +1 to the number of nodes and the number of arcs if your node ids start with 1.

Alternatively you can pass the proper number of nodes and arcs at the graph initialization and substract -1 from all of your node ids so they are zero based. But be careful – any arc is made of the node ids so you’ll have to substract -1 at the arcs level also.

Be the first to like.

April 10th, 2013

Solving the Minimum Cost Flow problem (5) – NetworkX

Network Simplex Solver: NetworkX

We will solve the instance of a Minimum cost flow problem described in (1) with NetworkX. Node 1 is the source node, nodes 2 and 3 are the transshipment nodes and node 4 is the sink node.

mcf_dimacs_1

Solving a Minimum Cost Flow Problem with NetworkX is pretty straight forward.

# import networkx
import networkx as nx

# create directed graph
G = nx.DiGraph()

# add node to graph with negative (!) supply for each supply node 
G.add_node(1, demand = -4)

# you can ignore transshipment nodes with zero supply when you are working with the mcfp-solver 
# of networkx
# add node to graph with positive (!) demand for each demand node
G.add_node(4, demand = 4)

# add arcs to the graph: fromNode,toNode,capacity,cost (=weight)
G.add_edge(1, 2, capacity = 4, weight = 2)
G.add_edge(1, 3, capacity = 2, weight = 2)
G.add_edge(2, 3, capacity = 2, weight = 1)
G.add_edge(2, 4, capacity = 3, weight = 3)
G.add_edge(3, 4, capacity = 5, weight = 1)

# solve the min cost flow problem
# flowDict contains the optimized flow
# flowCost contains the total minimized optimum
flowCost, flowDict = nx.network_simplex(G)

print "Optimum: %s" %flowCost  #should be 14

As you can see the code is very readable and easy to understand. The only thing you could get trouble with in the formulation of a min cost flow problem with networkx is the fact that the supply nodes get a negative supply value (because the python attribute is called „demand“) while the demand nodes require positive demand values.
This differs from the usual mathematical notation of supply and demand values in the network flow problems where supply values are positive and demand values of negative sign.
Also note that you only need to add supply and demand nodes to the graph. You don’t have to care of the transshipment nodes.

5 Leute mögen diesen Artikel.

April 6th, 2013

Solving the Minimum Cost Flow problem (4) – PuLP

Linear program solvers: PuLP

We will solve the instance of a Minimum cost flow problem described in (1) now with another linear program solver: PuLP. Node 1 is the source node, nodes 2 and 3 are the transshipment nodes and node 4 is the sink node.

mcf_dimacs_1

While lpsolve has this nice feature of reading DIMACS network flow problems, PuLP has nothing comparable to offer. So we have to transform the whole network flow problem into a plain linear program on our own.

At the PuLP wiki I found a sample for solving the transshipment problem with PuLP. Actually these code snippets work also for the minimum cost flow problem but in the original code there is a typo:

Wrong:

for n in Nodes:
    prob += (supply[n]+ lpSum([route_vars[(i,j)] for (i,j) in Arcs if j == n]) >=
             demand[n]+ lpSum([route_vars[(i,j)] for (i,j) in Arcs if i == n])), \
            "Steel Flow Conservation in Node:%s"%n

Right:

for n in Nodes:
    prob += (supply[n]+ lpSum([route_vars[(i,j)] for (i,j) in Arcs if j == n]) >=
             demand[n]+ lpSum([route_vars[(i,j)] for (i,j) in Arcs if i == n])), \
            "Steel Flow Conservation in Node %s"%n

The problem in the original code is the colon (‚:‘) at 'Steel Flow Conservation in Node:%s"%n'. This causes the generation of an incorrect *.lp file that is generated when you use PuLP. This *.lp file will be sent to a solver package of your (or PuLP’s) choice. So just remove the colon in this line.

Here we go with the full working code for our sample instance of the MCFP:

'''
Minimum Cost Flow Problem Solver with PuLP
Adapted from:
https://code.google.com/p/pulp-or/wiki/ATransshipmentProblem
The American Steel Problem for the PuLP Modeller
Authors: Antony Phillips, Dr Stuart Mitchell   2007

'''

# Import PuLP modeller functions
from pulp import *

# list of nodes
nodes = [1,2,3,4]

# supply or demand of nodes
            #NodeID : [Supply,Demand]
nodeData = {1:[4,0],
            2:[0,0],
            3:[0,0],
            4:[0,4]}

# arcs list
arcs = [ (1,2),
         (1,3),
         (2,3),
         (2,4),
         (3,4)]

# arcs cost, lower bound and capacity
            #Arc : [Cost,MinFlow,MaxFlow]
arcData = { (1,2):[2,0,4],
            (1,3):[2,0,2],
            (2,3):[1,0,2],
            (2,4):[3,0,3],
            (3,4):[1,0,5] }

# Splits the dictionaries to be more understandable
(supply, demand) = splitDict(nodeData)
(costs, mins, maxs) = splitDict(arcData)

# Creates the boundless Variables as Integers
vars = LpVariable.dicts("Route",arcs,None,None,LpInteger)

# Creates the upper and lower bounds on the variables
for a in arcs:
    vars[a].bounds(mins[a], maxs[a])

# Creates the 'prob' variable to contain the problem data    
prob = LpProblem("Minimum Cost Flow Problem Sample",LpMinimize)

# Creates the objective function
prob += lpSum([vars[a]* costs[a] for a in arcs]), "Total Cost of Transport"

# Creates all problem constraints - this ensures the amount going into each node is 
# at least equal to the amount leaving
for n in nodes:
    prob += (supply[n]+ lpSum([vars[(i,j)] for (i,j) in arcs if j == n]) >=
             demand[n]+ lpSum([vars[(i,j)] for (i,j) in arcs if i == n])), \
            "Flow Conservation in Node %s"%n

# The problem data is written to an .lp file
prob.writeLP("simple_MCFP.lp")

# The problem is solved using PuLP's choice of Solver
prob.solve()

# The optimised objective function value is printed to the screen    
print "Total Cost of Transportation = ", value(prob.objective)

The code is well commented. For a very detailed description have a look at the wiki of PuLP. Because there is no special network component in the PuLP modeler – any problem to solve is a linear program. So you have to write any instance of a MCFP as a plain LP which is not easy to understand at first glance.

3 Leute mögen diesen Artikel.

April 6th, 2013

Solving the Minimum Cost Flow problem (3) – lpsolve

Linear program solvers: Lpsolve

In the following articles we will solve the instance of a Minimum cost flow problem described in (1). Node 1 is the source node, nodes 2 and 3 are the transshipment nodes and node 4 is the sink node.

mcf_dimacs_1

As lpsolve can read the DIMACS file format we don’t have to transform a graph structure into a plain linear program. The XLI (eXtensible language interface) does this for us.
A good description of the DIMACS Minimum cost flow problem file format can be found here.
So for solving our MCFP we just have to write down the MCFP in DIMACS file format with a text editor of our choice (I called the file ’simple_MCFP.net‘) and let lpsolve do the dirty work:


c Problem line (nodes, links)
p min 4 5
c
c Node descriptor lines (supply+ or demand-)
n 1 4
n 4 -4
c
c Arc descriptor lines (from, to, minflow, maxflow, cost)
a 1 2 0 4 2
a 1 3 0 2 2
a 2 3 0 2 1
a 2 4 0 3 3
a 3 4 0 5 1

Here is the python code for solving a minimum cost flow problem in DIMACS file format with the python API of lpsolve:

from lpsolve55 import *

# full file path to DIMACS file
path = 'simple_MCFP.net'

# read DIMACS file and create linear program
lp = lpsolve('read_XLI', 'xli_DIMACS', path)

# solve linear program
lpsolve('solve', lp)

# get total optimum
optimum = lpsolve('get_objective', lp)

print optimum #should be 14

When you are working with the python API of lpsolve you’ll always have to call the function lpsolve() and pass the C method you want to access as first argument as string, e.g. lpsolve(’solve‘, lp). So obviously this API is a very thin layer upon the C library underneath.

The usage of the python API of lpsolve is not really pythonic because lpsolve is written in C – but hey: it works!  The complete list of functions can be found here.

Be the first to like.

März 31st, 2013

Solving the Minimum Cost Flow Problem (2) – Software

I don’t want to give you a complete overview of MCFP solvers because I just dipped into the world of linear and network flow programming. Instead of this I just want to mention a few MCFP solvers that I have tested. Python will be the programming language for the test – so I checked solvers with the corresponding language bindings.

Linear program solvers

As mentioned in the previous article one can solve MCFP also as ordinary linear programs. So here are the two packages I tested.

 

lpsolve is a Mixed Integer Linear Programming (MILP) solver written in C released under the LGPL. You can use lpsolve on Windows and Linux systems. For windows there is a nice GUI which is called LPSolveIDE. It reads several LP file formats such as *.mpl, *.lp and also the DIMACS file format which is used for network flow problems. Furthermore it has a rich API: you can use lpsolve with .NET (there is even a lpsolve plugin for Microsofts Solver Foundation MSF!), Java, PHP and Python.

 

PuLP is kind of a python wrapper for various (LP) solver packages released under the MIT license:

 PuLP is an LP modeler written in python. PuLP can generate MPS or LP files and call GLPK (http://www.gnu.org/software/glpk/glpk.html), COIN (http://www.coin-or.org/), CPLEX (http://www.cplex.com/), and GUROBI (http://www.gurobi.com/) to solve linear problems.

 

Specialized MCFP solvers

As the MCFP is a special linear program that can be solved more efficiently with algorithms that exploit the network structure you should also have a look at specialized solvers around.

 

NetworkX is a pure Python library that deals with graphs:

NetworkX is a Python language software package for the creation, manipulation, and study of the structure, dynamics, and functions of complex networks (http://networkx.github.com/).

The API is not very hard to learn, it’s well documented and obviously there is a active user community. NetworkX is released under the BSD license. It offers algorithms for many graph related issues (complete list here):

  • Clustering
  • Shortest Paths (Dijkstra, A*, Bellman-Ford, Floyd-Warshall)
  • Network flows (Minimum cost flow, maximum cost flow, minimum cut, Maximum flow minimum cost)
  • Minimum spanning tree

 

Google or-tools are a set of tools that deal not only with graph structures and algorithms but also with various other issues related to Operations Research (OR). This package is released under the Apache 2.0 license:

This project hosts operations research tools developed at Google and made available as open source. The project contains several tools:

  • A Constraint Programming solver.
  • A wrapper around third-party linear solvers (GLPK, CLP, CBC, SCIP, Sulum).
  • Knapsack algorithms.
  • Graph algorithms (shortest paths, min-cost flow, max flow, linear sum assignment).

Everything is coded in C++ and is available through SWIG in Python, Java, and .NET (using Mono on non-Windows platforms) (http://code.google.com/p/or-tools/).

Ok that’s it for now – let’s dive into the API of the packages in the next articles..

1 andere Person mag diesen Artikel.

März 24th, 2013

Solving the Minimum Cost Flow Problem (1) – Intro

Minimum cost flow problem: a short description

The Minimum cost flow problem (MCFP) is a generalized network flow problem that describes the goal of minimizing the total cost of flow through a network made of arcs and nodes which contains source, sink and transshipment nodes.

The minimum-cost flow problem is to find the cheapest possible way of sending a certain amount of flow through a flow network. Solving this problem is useful for real-life situations involving networks with costs associated (e.g. telecommunications networks), as well as in other situations where the analogy is not so obvious, such as where to locate warehouses (http://en.wikipedia.org/wiki/Minimum-cost_flow_problem).

Source nodes emit flow while sink nodes absorb flow. The flow through transshipment nodes is only intermediate: inflow must be equal to outflow at these nodes: so the net flow at these nodes is zero. Additionally there are maximum flow capacities on each arc of the network, i.e. it is not allowed that the amount of flow through a certain arc exceeds the maximum capacity specified on this arc.

In the following picture node 1 has 4 units of supply (= source node with positive value), node 4 has 4 units of demand (sink node with negative value) while nodes 2 and 3 are only transshipment nodes with a supply of zero. On the arcs you find lower flow bounds, upper flow bounds (or capacities) and costs. The goal is to ship 4 units from node 1 to node 4 with the lowest costs and without violation of the flow constraints (lower and upper bounds on the arcs).

Instance of a Minimum Cost Flow Problem

Instance of a Minimum Cost Flow Problem (http://lpsolve.sourceforge.net/5.5/DIMACS_mcf.htm)

The Minimum cost flow problem can be seen as a general problem from which various other network flow problems can be derived:

  • Transportation problem
  • Transshipment problem
  • Assignment problem
  • Shortest path (tree) problem

A special case of the Minimum Cost Flow Problem is the famous transportation problem (TP). In a regular TP the network lacks transshipment nodes and capacities on the arcs. So it’s a simpler version of the MCFP. The transshipment problem comes with transshipment nodes but lacks capacities on the arcs of the network. The assignment problem is a special version of the transportation problem and you can even solve a shortest path (tree) problem as a MCFP. Read more on the theoretical background here.

So in short: if you have a MCFP-Solver at hand you are able to solve various problems with this package.

The assignment problem is a special case of the transportation problem, which is a special case of the minimum cost flow problem, which in turn is a special case of a linear program. While it is possible to solve any of these problems using the simplex algorithm, each specialization has more efficient algorithms designed to take advantage of its special structure (http://en.wikipedia.org/wiki/Assignment_problem).

As stated in this quote and in previous articles, you can describe and solve a MCFP as a Linear program with an ordinary LP solver package. The simplex algorithm will likely be applied in this case. Or you can speed up things (magnitudes faster!) with a solver that is specialized on network flow problems and exploits the network structure of the MCFP: a network simplex solver or the push-relabel algorithm.

More on this in the next articles..

2 Leute mögen diesen Artikel.

Januar 20th, 2012

Geoprocessing mit Open Source GIS Tools und Python

Aufgabenstellung

Momentan soll ich eine ganze Menge an Geodaten für einen wissenschaftlichen Auftraggeber verarbeiten. Dabei geht es um die Potentialabschätzung der Wuchsleistung von schnellwachsenden Baumarten. Ich freute mich riesig auf den professionellen Einsatz von Open Source GIS Tools unter etwas anspruchsvolleren Bedingungen – denn  in letzter Zeit hatte ich beruflich wenig mit dem OS GIS Universum zu tun.

Große Datensätze

Anspruchsvoll? Oh ja – es handelt sich um eine bayernweite Analyse und da kommt schnell einiges zusammen…
Hier eine kleine Übersicht:

  • 4 Sachdatenbanken (MS Access)
  • 9 vektorielle Geodatensätze in insgesamt über 750 einzelnen ESRI Shapefiles (ca. 6,8 GB)
  • 3 Rasterdatensätze als ESRI ASCII bzw. Binary GRID (ca. 4,6 GB)

Nicht übel. Diese Menge von rund 11 GB sollte also mit Open Source GIS Tools verarbeitet werden. Das war im Auftrag natürlich nicht gefordert, da dieser softwareunabhängig formuliert war – aber für mich war es klar, dass ich FOSSGIS-Tools einsetze.

Ehrlich gesagt habe ich bislang mit FOSSGIS-Werkzeugen keine aufwendigen GIS-Analysen durchgeführt, sondern eher ausprobiert, was damit alles geht. Nun also eine vollständige Analyse. Dass das spannend wird, dachte ich mir schon…

Welche Werkzeuge?

Nun gut – welches Werkzeug solls denn nun sein? Da die Daten – vor allem die Shapefiles – ziemlich zerhackstückelt über ganz Bayern daher kamen (Schönen Gruß ans LfU und das LVG), mussten diese erstmal zu bayernweiten Datensätzen zusammengeführt werden. Das war vor Angebotsabgabe leider nicht bekannt – aber gut, ich nahm die Herausforderung an.

QGIS – fTools

Bei QGIS gibts eine Erweiterung, welche auch Shapefiles zusammenführen kann. Die sog. fTools können dies und auch noch viel mehr (z.B. Intersection, Dissolve, Multipart to single part etc.). Dabei basiert fTools auf QGIS und GDAL/OGR, um die Geoprocessing-Funktionen anbieten zu können. Außerdem ist es in Python geschrieben.
Die Variante der fTools in QGIS kam aufgrund der doch langen Laufzeiten nicht in Frage. Die GUI von QGIS frisst da anscheinend doch einiges an Performance — außerdem muss fTools dann als Plugin von QGIS mit dem Mutterprozess hin und her kommunizieren. Das wollt ich mir ersparen.
Man muss aber sagen, dass die fTools für kleinere Batch-Arbeiten gut geignet sind. Das Pythonskript, das ich mit der Vorlage der fTools geschrieben habe liest sich einfach. Das QGIS-Objektmodell folgt dem üblichen „Feature – Geometry – Attributtable“-Schema und ist somit leicht nachzuvollziehen.
Bis die QGIS-Pythonbindings dann aber mal ansprechbar sind kann ne Zeit vergehen. Gary Sherman hat dazu eine ausführliche Anleitung hier. Dazu mal an anderer Stelle mehr.

GDAL/OGR – kennste eine – kennste alle?

GDAL/OGR ist ja das ETL (Extract-Transform-Load) Werkzeug schlechthin. Sagt man. Ist auch so sag ich jetzt, aber noch vor ein paar Wochen hatte ich noch nicht viel Ahnung, wie man die Bibliothek einsetzt.
Ich hatte GDAL/OGR schon immer unbedacht benutzt, wenn ich z.B. mit QGIS oder GRASS arbeitete, da viele Desktop-GIS diese Bibliothek unter der Haube haben – allerdings eben immer „nur“ als Datenprovider für Desktop GIS.

Nun war klar, dass ich die Python-Bindings von GDAL/OGR genauer anschauen muss. Also erstmal recherchieren und viel ausprobieren (die Python-Doku für GDAL/OGR ist derzeit noch sehr spärlich) und irgendwann klappts dann aber auch. Das grundsätzliche Modell bei Geoprocessing-Bibliotheken ist ja immer relativ ähnlich: es gibt ein Feature-Objekt mit Sachdaten und ner Geometrie und man hat ein paar GIS-Methoden (welche im FOSSGIS-Universum ja eigentlich fast immer auf GEOS und somit der JTS basieren) zur Verfügung (Intersection, Buffer, etc.). Wie soll es auch mordsmäßig verschieden sein?
Genau so ists auch bei den Python-Bindings von GDAL/OGR. Die Extract und Load-Komponente kommt von GDAL/OGR selbst. Die geometrischen Verschneidungsoperatoren werden bei GDAL/OGR mit GEOS bewerkstelligt, während Transformationen mit PROJ4 behandelt werden.

Wer – wie ich – eher aus der ESRI-Ecke mit dem ArcGIS Geoprocessor kommt, wird sich bei GDAL/OGR in Python auch recht schnell zurecht finden. Die FDO API in C# war da schon ne andere (weil abstraktere) Nuss.

Interessante Links für GDAL/OGR mit Python:

Letztlich konnten die Daten dann mit einem  Pythonskript und der GDAL/OGR-Bibliothek zu bayernweiten Datensätzen zusammengeführt werden. Und wie schnell das geht! Da ist GDAL/OGR wirklich in ihrem Element. Wahnsinn, was die Bibliothek wegschreibt:

  • rund 100 Shapefiles mit insgesamt 2,1 Mio Features in 20 Minuten (800MB)
  • rund 100 Shapefiles mit insgesamt 3,1 Mio Features in 35 Minuten (4,2 GB)

Diese Zahlen überzeugen.

PostGIS

Auch PostGIS fällt nicht aus der Reihe der FOSSGIS-Tools und basiert auf GEOS und der PROJ4-Bibliothek. Die Kombination von SQL als Datenabfragesprache mit geometrischen Funktionen (Spatial SQL) fand ich immer spannend. Nachdem die Shapefiles also zusammengeführt waren, wollte ich einfach alle Vektordaten in eine PostGIS-DB importieren und die Aufbereitung in PostGIS machen. Auch dass mit PL/Python meine Lieblingssprache Python in einer Datenbank als prozedurale Sprache zur Verfügung steht find ich klasse (ich denke hier an die Lesbarkeit von PL/SQL Code vs. Python…).
Allerdings kam ich nicht weit. PostGIS hat definitiv als professionelle Geodatenhaltung seine Stärken und einige Anfragen an Daten kann es auch schnell beantworten, aber PostGIS ist nicht „post GIS“ (im Sinne von „nach GIS“) wie es im Vorwort des Buchs „PostGIS in Action“ hies. Man braucht bei solchen Datenmengen anscheinend schon noch ein extra Analysewerkzeug. Vielleicht liegts auch an meiner Rechenpower oder meinen PostGIS/SQL-Kenntnissen, aber GIS-Analysen mit dieser Datenmenge klappt nach meiner Erfahrung nicht in PostGIS.  Nach Laufzeiten von 2 Tagen habe ich die aufwendigen Verschneidungen (bayernweit) sein gelassen. Schade eigentlich. Ich hätte irgendwie mehr Performance von PostGIS erwartet.

SAGA GIS

Raster? Nimm SAGA GIS. Eben – für Rasterdaten ist SAGA GIS spitze. Performant, gut bedienbar, sogar per Python oder Kommandozeile ansprechbar für Batchabläufe. Bei Vektordaten hat SAGA aber so seine Macken. Außerdem werden nicht viele Möglichkeiten in der Vektoranalyse angeboten. Diejenigen, die dabei sind basieren interessanterweise jedoch nicht auf GEOS, sondern auf CGAL.
Aber gut, SAGA GIS ist eben vor allem ein Raster-GIS. In Ordnung.

GRASS GIS – hat so nen Bart?

Allerdings. GRASS GIS hat schon einige Jahre auf dem Buckel (30+?) – wird aber laufend fortentwickelt. Seit geraumer Zeit sieht auch die GUI etwas ansprechender aus. Von den Analysemöglichkeiten her ist GRASS GIS im FOSSGIS-Universum meiner Meinung nach einsame spitze.

Nimm GRASS GIS – wenn Du nicht ewig Zeit hast. Im Zusammenhang mit Scripting ist das Werkzeug ein verdammt performantes Analyse-Tool. Natürlich sollte man an dieser Stelle die etwas höhere Lernkurve bei GRASS ggü. anderen Werkzeugen berücksichtigen. Aber die Zeitinvestition lohnt sich!

Ich bin in den o.g. Ausführungen immer auf die Geometriebibliothek der Tools eingegangen. Da fast alle OS GIS Werkzeuge letztlich die geometrischen Funktionen an die GEOS-Bibliothek weiterreichen, ist auch die Performance bei großen Datenmengen vergleichbar.
Damit ist es relativ egal, ob man GDAL/OGR oder die fTools von QGIS oder PostGIS verwendet. Alle basieren auf GEOS. GEOS ist ja ein C++ Port der JTS (Java Topology Suite). Und diese ist wiederum in anderen FOSSGIS-Tools verbaut wie gvSIG, OpenJUMP, uDIG oder Kosmo. GIS-Analysen mit derartigen Datenmengen in Java-Tools durchzuführen hatte bei mir leider regelmäßig einen „Out of heap space“ Fehler zur Folge. Deshalb schied auch die JTS aus.

Die GRASS GIS Overlay-Operatoren sind jedoch als C-Algorithmen anscheinend richtig speicheroptimiert geschrieben. Trotzdem bleibt ein kleiner Wermutstropfen: Auch in GRASS GIS konnte ich die Datenmengen nicht bayernweit importieren bzw. verschneiden. Mit einem 64bit System und mehr RAM wäre das vielleicht schon gegangen. So aber blieb mir nur ein „Tiled Processing“, in dem ich die bayernweiten Datensätze in 16 Kacheln aufgeteilt habe. Mit GRASS GIS hats dann im Batch-Modus per Pythonbindings allerdings recht flott (rund 1,5 Std. für einen bayernweiten Verschneidungsdurchlauf) geklappt. Vielleicht mach ich demnächst auch nochmal nen Vergleich mit einem Tiled Processing der Daten in PostGIS, aber ich habe wenig Hoffnung, dass dies ähnlich performant ist wie GRASS GIS.

Und ne kleine GRASS Python Utility Klasse ist auch noch herausgekommen:

# -*- coding: latin-1 -*-
#  Helper for GRASS GIS Python Scripting
#  Johannes Sommer, 18.01.2012

class utility():
	import core as grass
	''' Utility-Class for GRASS GIS Python Scripting.
		Copyright by Johannes Sommer 2012.

		- doAddColumn
		- doCalcColumnValue
		- doRenameColumn
		- doIntersection
		- doDifference
		- getMapsets
		- getVectorLayers
		- getRasterLayers
		- getVectorLayersFromMapset
		- getRasterLayersFromMapset
		- getColumnList
		- doFilterLayerList 

		Usage:
		import grass.utils as utils
		gUtil = util.utility()
		gUtil.dodifference('forest', lakes')
		'''

	def doAddColumn(self, _lyrA, _colName, _colType):
		''' Adds a column to the input vector layer with the given type.

			GRASS datatypes depend on underlying database, but all
			support:
			- varchar (length)
			- double precision
			- integer
			- date

			Usage:
			doAddColumn('lakes', 'width', 'double precision')	'''
		retVal = grass.run_command('v.db.addcol', map=_lyrA,
								   columns='%s %s' % (_colName, _colType) )
		return retVal

	def doCalcColumnValue(self, _lyrA, _colName, _newVal):
		''' Updates a column of the input vector layer with the passed value.
			@_newVal can be a certain value or a calculation of existing columns.

			Usage:
			doCalcColumnValue('lakes', 'width', 'length/2')		'''
		retVal = grass.run_command('v.db.update', map=_lyrA,
								   column=_colName, value=_newVal)
		return retVal

	def doRenameColumn(self, _lyrA, _colNameA, _colNameB):
		''' Renames a column of the input vector layer.

			Usage:
			doRenameColumn('lakes', 'width', 'width_old')		'''
		retVal = grass.run_command('v.db.renamecol', map=_lyrA,
								   column='%s,%s' %(_colNameA,_colNameB))
		return retVal

	def doIntersection(self, _lyrA, _lyrB, OVERWRITE_OUTPUT==False):
		''' Performs a geometric intersection ('AND'-Operator) of two vector layers
			and writes the output to a layer called 'layerNameA_X_layerNameB'.

			@OVERWRITE_OUTPUT can be set to 'Y' or 'N'
			Usage:
			doIntersection('lakes', 'forest')	'''
		outputLyr = _lyrA.split('@')[0] + '_X_' + _lyrB.split('@')[0]

		optArg = ''
		if OVERWRITE_OUTPUT == True:
			optArg = '--overwrite'
		retVal = grass.run_command('v.overlay %s' % optArg, ainput=_lyrA,
								   binput=_lyrB, output=outputLyr,
								   operator='and')
		return retVal

	def doDifference(self, _lyrA, _lyrB, OVERWRITE_OUTPUT==False):
		''' Performs a geometric difference ('NOT'-Operator) of two vector layers
			and writes the output to a layer called 'layerNameA_DIFF_layerNameB'.

			@OVERWRITE_OUTPUT can be set to 'Y' or 'N'
			Usage:
			doDifference('lakes', 'forest')	'''
		outputLyr = _lyrA.split('@')[0] + '_DIFF_' + _lyrB.split('@')[0]

		optArg = ''
		if OVERWRITE_OUTPUT == True':
			optArg = '--overwrite'
		retVal = grass.run_command('v.overlay %s' % optArg, ainput=_lyrA,
								   binput=_lyrB, output=outputLyr,
								   operator='not')
		return retVal

	def getMapsets(self):
		''' Returns all GRASS mapsets in current GRASS location as list.'''
		return grass.mapsets()

	def getVectorLayers(self):
		''' Returns all vector layers in current GRASS location as list.'''
		return grass.list_strings('vect')

	def getRasterLayers(self):
		''' Returns all raster layers in current GRASS location as list.'''
		return grass.list_strings('rast')

	def getVectorLayersFromMapset(self, _mapset):
		''' Returns all vector layers in given mapset as list.'''
		vLyrList = getVectorLayers()
		# Filter List by MAPSET
		vLyrListFiltered = []
		for item in vLyrList:
			if _mapset in item:
				vLyrListFiltered.append(item)
		return vLyrListFiltered

	def getRasterLayersFromMapset(self, _mapset):
		''' Returns all raster layers in given mapset as list.'''
		rLyrList = getRasterLayers()
		# Filter List by MAPSET
		rLyrListFiltered = []
		for item in rLyrList:
			if _mapset in item:
				rLyrListFiltered.append(item)
		return rLyrListFiltered

	def getColumnList(self, _vLyr):
		''' Returns all column names of a vector layer as list.'''
		# Import vector wrapper from %GISBASE%\etc\python\script
		import vector as vector
		out_colList = []
		out_colList = vector.vector_columns(_vLyr, getDict = False)
		return out_colList

	def doFilterLayerList(self, _lyrList, _lyrkey):
		''' Filters input layer list per _lyrkey and returns filtered list.'''
		out_lyrList = []
		for item in _lyrList:
			if _lyrkey in item:
				out_lyrList.append(item)
		return out_lyrList

9 Leute mögen diesen Artikel.

Januar 12th, 2011

ArcGIS Server REST-Schnittstelle mit Python auslesen

Wenn mit ESRI ArcGIS Server Dienste veröffentlicht werden, kann für jeden Dienst die REST-Schnittstelle abgerufen werden.
Diese offenbart einige Informationen über den angebotenen Dienst. Sie kann ganz simpel mit dem Browser angesteuert werden wie die Dokumentation von ESRI verrät:

The ArcGIS Server REST API, short for Representational State Transfer, provides a simple, open Web interface to services hosted by ArcGIS Server. All resources and operations exposed by the REST API are accessible through a hierarchy of endpoints or Uniform Resource Locators (URLs) for each GIS service published with ArcGIS Server.

When using the REST API, you typically start from a well-known endpoint, which represents the server catalog. The table of contents for this Help system mimics the hierarchy of resources in the REST API. The default start URL for an ArcGIS Server installation is:

  • Java:

http://<host>:8399/argis/rest

  • .NET:

http://<host>/arcgis/rest

Ein Beispiel für einen Aufruf der REST-Schnittstelle wäre demnach:

http://server.arcgisonline.com/arcgis/rest/services

Das Ergebnis im Browser entspricht dem Wurzelverzeichnis eines ArcGIS Servers. In dieser hierachischen Struktur kann man sich nun manuell durchklicken um zu Informationen über die angebotenen Dienste und ihrer Layer zu bekommen.

Dabei werden die Elementbezeichnung wie Folder oder Layernamen bei jedem Klick in die angezeigte URL übernommen – nach dem Schema:

http://<host>/arcgis/rest/services/<folder>/<service>/MapServer/

Ein beispielhafter Aufruf zu einem Layer im Dienst USA_1990-2000_Population_Change wäre damit:

http://server.arcgisonline.com/ArcGIS/rest/services/Demographics/USA_1990-2000_Population_Change/MapServer/

Hängt man an den Aufruf eine Formatanweisung hinzu, kann die Information nicht nur hübsch aufbereitet als HTML-Seite bezogen werden, sondern auch als JSON-Objekt.
Der Schlüssel dazu lautet: ?f=pjson und die gesamte URL nun:

http://server.arcgisonline.com/ArcGIS/rest/services/Demographics/USA_1990-2000_Population_Change/MapServer/?f=pjson

Ein Ausschnitt des Ergebnisses diesen Aufrufs sieht folgendermaßen aus:

{
  "currentVersion" : 10.01,
  "serviceDescription" : "This thematic map indicates the annual compound rate of total population change in the United States from 1990 to 2000. Total Population is the total number of residents in an area.  Residence refers to the “usual place” where a person lives.  Total Population for 2000 is from the U.S. Census 2000.  Total Population for 1990 is from the 1990 U.S. Census adjusted to 2000 boundaries.  The geography depicts states at greater than 25m scale, counties at 1m to 25m scale, Census Tracts at 250k to 1m scale, and Census Block Groups at less than 250k scale. The map has been designed to be displayed with semi-transparency of about 50% for overlay on other base maps, which is reflected in the legend for the map.  For more information on this map, visit us \u003ca href=\"http://goto.arcgisonline.com/maps/Demographics/USA_1990-2000_Population_Change \" target=\"_new\"\u003eonline\u003c/a\u003e.",
  "mapName" : "Layers",
  "description" : "This thematic map indicates the annual compound rate of total population change in the United States from 1990 to 2000. Total Population is the total number of residents in an area. Residence refers to the “usual place” where a person lives. Total Population for 2000 is from the U.S. Census 2000. Total Population for 1990 is from the 1990 U.S. Census adjusted to 2000 boundaries. The geography depicts states at greater than 25m scale, counties at 1m to 25m scale, Census Tracts at 250k to 1m scale, and Census Block Groups at less than 250k scale. The map has been designed to be displayed with semi-transparency of about 50% for overlay on other base maps, which is reflected in the legend for the map.  For more information on this map, visit us online at http://goto.arcgisonline.com/maps/Demographics/USA_1990-2000_Population_Change",
  "copyrightText" : "Copyright:© 2010 ESRI",
  "layers" : [
    {
      "id" : 0,
      "name" : "USA 1990-2000 Population Change",
      "parentLayerId" : -1,
      "defaultVisibility" : true,
      "subLayerIds" : [1, 2, 3, 4],
      "minScale" : 0,
      "maxScale" : 0
    },
    {
      "id" : 1,
      "name" : "Block Groups",
      "parentLayerId" : 0,
      "defaultVisibility" : true,
      "subLayerIds" : null,
      "minScale" : 150000,
      "maxScale" : 0
     }

Dieses JSON-Objekt kann man nun in eigenem Code, der die URL ausliest, wunderbar weiterverarbeiten. Wunderbar deswegen, weil eine HTML-Ausgabe auszulesen dagegen ein ziemliches Gefrickel ist.

Ich wollte das Ganze mit Python anzapfen um Layerinformationen aus dem ArcGIS Server Dienst zu extrahieren. Dazu braucht es neben der urllib2 von Python noch das Modul simplejson, um die Informationen, die im JSON-Objekt gekapselt sind, auf einfache Weise in Python zugänglich zu machen.

Nach der Installation von simplejson kann es losgehen. Wir importieren simplejson und das Standardmodul urllib2, das die Verbindung zum Server übernimmt.

import simplejson as json
import urllib2
URL = "http://server.arcgisonline.com/ArcGIS/rest/services/Demographics/"

Anschließend definiere ich eine Funktion, die jeweils die Layer-ID und den Layernamen aus dem angegebenen Dienst anzeigen soll. Die Funktion benötigt nur eine URL zu einem ArcGIS Server und den Servicenamen, der abgefragt werden soll.  Selbstverständlich könnte man auch ganz oben in der Hierachie beginnen und aus dem Wurzelverzeichnis alle Diensteverzeichnisse auflisten lassen und dann weiter in die jeweiligen Dienste schauen. An dieser Stelle möchte ich jedoch nur einen einfachen Fall beschreiben.

def getLayerInfoFromREST(url, _serviceName):
    ''' Liest die Layer-IDs des ArcGIS REST Dienstes aus
        der angegebenen URL und liefert eine Liste mit den
        Layer-IDs (Strings) und den Layernamen zurück.
        Benötigtes Format der URL:
        http://<host>/ArcGIS/rest/services/<folder>/'''

    serviceName = _serviceName
    url = url + serviceName + '/MapServer/?f=pjson'

    response = urllib2.urlopen(url).read()
    jsonObj = json.loads(response)

Wir haben also die URL mit dem Servicenamen angereichert, eine Antwort des Servers auf unseren Request in das Objekt response gespeichert und zu guter Letzt die response als JSON-Objekt geladen. Denn auf unsere Anfrage in der URL mit der Endung ?f=pjson bekommen wir auf jeden Fall ein JSON-Objekt zurück. Jetzt initialisieren wir das Dictionary resultDict, das unser Ergebnis der Abfrage speichern wird. Auch das Element jsonObj ist ein Dictionary.
Den Inhalt des JSON-Objekts durchlaufen wir mit einer for-Schleife und der Funktion iterkeys(). Genauer gesagt durchlaufen wir nur die Keys des Dictionarys. Wenn der Key bzw. anschaulicher der Knoten „layers“ gefunden wurde, steigen wir mit der zweiten for-Schleife in diesen Knoten ein und iterieren über die untergeordneten Elemente. Auf dieser Ebene können nun die Layerinformationen abgerufen werden, z.B. die ID und der Name, welche mit der Funktion get() und dem entsprechendem Key abfragbar sind.
Nun wird das Ergebnis in das Dictionary resultDict gesichert.

    resultDict = {}
    for key in jsonObj.iterkeys():
        # layers entspricht dem Data Frame Namen des zugrunde-
        # liegenden ArcMap-Dokument und muss evtl. angepasst
        # werden
        if key == 'layers':
            root = jsonObj.get(key)

            for layerInfo in root:
                nameString = layerInfo.get('name')
                id = int(layerInfo.get('id'))

                resultDict[id] = nameString

    return resultDict

Jede Funktion muss auch aufgerufen werden, und genau das machen wir im folgenden. Die Parameter sind wie o.g. URL und ein Servicename.

# Funktionsaufruf
resultDict = getLayerInfoFromREST(URL, "USA_1990-2000_Population_Change")

for item in resultDict.iteritems():
    print item

In der Ausgabe des Skripts sehen wir nun die Layer-IDs und die Namen der Layer aus dem ArcGIS Server Dienst.

>>>
(0, 'USA 1990-2000 Population Change')
(1, 'Block Groups')
(2, 'Tracts')
(3, 'Counties')
(4, 'States')
>>>

Das Beispielsskript „jsonRESTSample.py“ kann hier gedownloaded werden.

5 Leute mögen diesen Artikel.

August 8th, 2010

Spaß mit Python und einem WFS (2)

Wie bereits erwähnt wollte ich den WFS Altimeter mit Python ansprechen und nicht über die Javamaske des BKG bzw. im Browser. Da man den WFS nur mit POST-Requests füttern kann ist letzteres ein mühsames Unterfangen. Deshalb musste der Dienst mit Python angesprochen werden.

Nach einer kurzen Suche fand ich OWSLib für Python. Es gibt also bereits ein Pythonmodul, dass sich dem Thema OGC Web Services (OWS) angenommen hat. Damit können WMS und WFS in Python angesprochen werden, ohne dass man sich durch die XMLs der Dienste hindurch parsen müsste.

weiterlesen »

1 andere Person mag diesen Artikel.