Kurs Python richtig lernen/Aufgabe 2: Raster Wohnungen - Versionsgeschichte

Stefan am 16. Januar 2013 um 12:43 Uhr

2013-01-16T12:43:52Z

← Nächstältere Version		Version vom 16. Januar 2013, 12:43 Uhr
Zeile 1:		Zeile 1:
			Zurück zum [[Kurs Python richtig lernen#Uebungen]].

	* Einlesen einer Textdatei mit Wohnungsinseraten (appartements.txt, enthalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]])		* Einlesen einer Textdatei mit Wohnungsinseraten (appartements.txt, enthalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]])
	* Koordinatentransformation von WGS 84 nach CH1903 (Schweizer Landeskoordinaten)		* Koordinatentransformation von WGS 84 nach CH1903 (Schweizer Landeskoordinaten)

Josh.x am 19. September 2012 um 14:00 Uhr

2012-09-19T14:00:44Z

← Nächstältere Version		Version vom 19. September 2012, 14:00 Uhr
Zeile 3:		Zeile 3:
	* Verwendung von Numpy-Arrays		* Verwendung von Numpy-Arrays
	* Schreiben einer Erdas-Imagine Datei (öffenbar in [[Quantum GIS\|QGIS]])		* Schreiben einer Erdas-Imagine Datei (öffenbar in [[Quantum GIS\|QGIS]])
			Installation benötigter Module in Linux (Terminal):
			sudo apt-get install python-gdal
			sudo apt-get install python-numpy
			sudo apt-get install python-pyproj

	== mietpreis.py ==		== mietpreis.py ==

Josh.x am 19. September 2012 um 13:47 Uhr

2012-09-19T13:47:33Z

← Nächstältere Version		Version vom 19. September 2012, 13:47 Uhr
Zeile 1:		Zeile 1:
			* Einlesen einer Textdatei mit Wohnungsinseraten (appartements.txt, enthalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]])
			* Koordinatentransformation von WGS 84 nach CH1903 (Schweizer Landeskoordinaten)
			* Verwendung von Numpy-Arrays
			* Schreiben einer Erdas-Imagine Datei (öffenbar in [[Quantum GIS\|QGIS]])

	== mietpreis.py ==		== mietpreis.py ==

	Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten ~~(appartements.txt, enthalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]])~~. Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.		Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten. Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.

	<nowiki>		<nowiki>

Josh.x: /* mietpreis.py */

2012-09-19T13:24:38Z

mietpreis.py

← Nächstältere Version		Version vom 19. September 2012, 13:24 Uhr
Zeile 1:		Zeile 1:
	== mietpreis.py ==		== mietpreis.py ==

	Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten (appartements.txt, ~~etnhalten~~ in [[Media:Musterloesungen.zip\|Musterloesungen.zip]]). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.		Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten (appartements.txt, enthalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]]). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.

	<nowiki>		<nowiki>

Josh.x: /* mietpreis.py */

2012-09-19T13:24:18Z

mietpreis.py

← Nächstältere Version		Version vom 19. September 2012, 13:24 Uhr
Zeile 1:		Zeile 1:
	== mietpreis.py ==		== mietpreis.py ==

	Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten (appartements.txt). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.		Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten (appartements.txt, etnhalten in [[Media:Musterloesungen.zip\|Musterloesungen.zip]]). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.

	<nowiki>		<nowiki>

Josh.x: Die Seite wurde neu angelegt: „== mietpreis.py == Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei…“

2012-09-19T12:52:48Z

Die Seite wurde neu angelegt: „== mietpreis.py == Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei…“

Neue Seite

== mietpreis.py ==

Berechnet den durchschnittlichen Mietpreis pro Quadratmeter auf einem Raster das über die ganze Schweiz geht. Ausgangspunkt ist die Textdatei mit den Wohnungsinseraten (appartements.txt). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.

<nowiki>
# -*- coding: utf-8 -*-

# Input-Parameter:
infile = 'appartements.txt'
outfile = 'appartements.img' # Raster-Datei im Erdas Imagine Format
env = 480000, 75000, 860000, 300000 # Envelope: xmin, ymin, xmax, ymax
res = 1000 # Resolution: Grösse einer Rasterzelle

import numpy as np
import geo # Dies ist unser eigenes Modul (geo.py)

# Raster-Grösse berechnen (wie viele Zellen in x und y?)
nx = (env[2] - env[0]) / res
ny = (env[3] - env[1]) / res

# Um den durchschnittlichen Mietpreis pro m2 zu berechnen brauchen wir folgende
# Informationen:
# surface_habitable (m2), loyer_mois (CHF)
# Wir müssen für jede Rasterzelle die Summe von CHF/m2 speichern, die Anzahl
# Wohnungen, und den Durchschnitt (Summe / Anzahl). Also brauchen wir drei
# Numpy-Arrays.
summe = np.zeros((nx,ny))
anzahl = np.zeros((nx,ny))
durchschnitt = np.zeros((nx,ny))

# Die Projektion vorbereiten.
wgs84_to_ch1903 = geo.Proj(from_srs=4326, to_srs=21781)

# Die Input-Datei öffnen
fp = open(infile)
h = fp.readline() # Den Header mit den Variablennamen lesen

# Zeile für Zeile der Input-Datei behandeln
for zeile in fp:
werte = zeile.split('\t')
# Falls wir weniger als 5 Werte haben, Zeile überspringen:
if len(werte) < 5: continue
# Falls eine Kolonne einen fehlenden Wert hat Zeile überspringen:
if werte[4] == 'NULL' or werte[5] == 'NULL': continue
lat, lng = float(werte[1]), float(werte[2])
flaeche = float(werte[4])
mietpreis = float(werte[5])
# Falls der Mietpreis über 20'000 Franken ist, haben wir es wahrscheinlich
# mit einer Falschinformation zu tun. Also Zeile überspringen
if mietpreis > 20000: continue
# Nun können wir die Koordinaten in Schweizer Koordinaten umwandeln:
x, y = wgs84_to_ch1903.transform(lng, lat)
# Und die dazugehörigen Pixel-Koordinaten finden:
px, py = geo.coord_to_pixel([x,y], env, res)
# Den Mietpreis pro m2 berechnen:
mietpreis_m2 = mietpreis / flaeche
# Und nun die drei Numpy-Arrays updaten:
summe[px, py] += mietpreis_m2
anzahl[px, py] += 1
durchschnitt[px, py] = summe[px, py] / anzahl[px, py]

# Wir sind durch die ganze Datei durch. Also Datei schliessen.
fp.close()

# Und nun den Numpy-Array in eine Erdas-Imagine Datei schreiben:
geo.write_raster(outfile, durchschnitt, env, res)
</nowiki>

== geo.py ==

Hilfsmodul für mietpreis.py.

<nowiki>
# -*- coding: utf-8 -*-

import pyproj
import numpy as np
import osgeo.gdal as gdal

class Proj:
"""
Eine Klasse die die Projektion von Koordinaten etwas vereinfacht.
Zuerst eine Instanz der Klasse erstellen mit den korrekten
Koordinatensystemen (EPSG-Code als String angeben), und dann die Koordinaten
projezieren:
wgs84_to_ch1903 = Proj(4326, 21781)
wgs84_to_ch1903.transform(6.1, 46.5)
"""
def __init__(self, from_srs, to_srs):
self.s_srs = pyproj.Proj(init='EPSG:'+str(from_srs))
self.t_srs = pyproj.Proj(init='EPSG:'+str(to_srs))

def transform(self, x, y):
return pyproj.transform(self.s_srs, self.t_srs, x, y)

def coord_to_pixel(coord, env, res):
"""
Transformiert die Geo-Koordinaten in Pixel-Koordinaten.
Dazu brauchen wir nicht nur die Koordinaten, sondern auch die Envelope
(Tupel mit xmin, ymin, xmax, ymax), und die Raster-Auflösung (res).
"""
x = int((coord[0] - env[0]) / res)
y = int((coord[1] - env[1]) / res)
return x, y

def write_raster(outfile, nparray, env, res):
"""
Schreibt einen Numpy-Array in eine Erdas-Imagine-Datei.
Die GeoTransform enthält eine Liste
[upper_left_x, resolution_x, rotation_1, upper_left_y, rotation_2, -resolution_y]
Die Rotations-Parameter sind in einem geraden Bild 0.
"""
driver = gdal.GetDriverByName('HFA')
nx, ny = nparray.shape
dataset = driver.Create(outfile, xsize=nx, ysize=ny,
bands=1, eType = gdal.GDT_Float32)
dataset.SetGeoTransform([env[0], res, 0.0, env[3], 0.0, -res])
band1 = dataset.GetRasterBand(1)
band1.WriteArray(np.flipud(nparray.T))
</nowiki>