Kurs Python richtig lernen/Aufgabe 2: Raster Wohnungen

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 Musterloesungen.zip). Dabei sind die Koordinaten der Wohnungen in Lat/Long, und wir konvertieren sie hier in Schweizer Koordinaten.

# -*- 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)
 

geo.py

Hilfsmodul für mietpreis.py.

# -*- 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))