Kurs Python richtig lernen/Aufgabe 2: Raster Wohnungen: Unterschied zwischen den Versionen
Aus Geoinformation HSR
Josh.x (Diskussion | Beiträge) (→mietpreis.py) |
Josh.x (Diskussion | Beiträge) (→mietpreis.py) |
||
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, | + | 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> |
Version vom 19. September 2012, 14:24 Uhr
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 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))