Gode Pythonpakker¶
I dette kapitel beskriver vi en række gode og anbefalelsesværdige Pythonpakker, der kan hjælpe med at løse gængse problemstillinger en udvikler i SDFI ofte vil støde på. I udvælgelsen af pakkerne herunder er der lagt vægt på særligt pythonske implementeringer, der benytter sig af koncepter introduceret i nyere udgaver af Python.
Håndtering af geodata¶
GDAL er det mest udbredte programmel til læsning, skrivning og manipulation af geodata. GDAL kan bearbejde både raster- og vektordata. Vektordata manipuleres med underbiblioteket OGR. GDAL-projektet tilbyder pythonpakke men den er ikke særligt intuitiv at bruge for pythonudviklere, da der er tale om en en-til-en mapning af det underliggende C API. Bedre alternativer er pakkerne Rasterio (rasterdata) og Fiona (vektordata), der begge tilbyder et bedre og mere pythonsk API til GDAL og OGR.
Rasterio¶
Rasterio bruges til at læse og skrive rasterdata. Herunder eksempler på begge dele.
Læsning af rasterdata er utroligt simpelt og minder meget om måden almindelige tekstfiler læses med Python:
import rasterio
with rasterio.open("eksempel.tif", "r") as grid:
# Læs information om hele filen
epsg_kode = grid.crs
gridstørrelse = grid.shape
# Behandl data for hvert enkelt bånd i rasterfilen.
# Bemærk at data returnes i form af numpy arrays.
for båndnr in grid.indexes:
data = grid.read(båndnr)
gennemsnit = data.mean()
Oprettelse af en ny fil og eftefølgende skrivning af nyt data er mere omfattende. Det forudsættes at der er taget stilling til filens dækningsområde, koordinatsystem, antal bånd, gridstørrelse, filformat og datatype. Alt dette kan dog klares relativt simpelt, som demonstreret i eksemplet herunder, hvor vi laver en tiff-fil med tre bånd, der omtrentligt dækker det danske område og registrerer data i ETRS89-koordinater.
import rasterio
import numpy as np
# Grid dimensioner
antal_bånd = 3
højde = 256
bredde = 256
# Dækningsområde
vest = 8.0
øst = 13.0
syd = 54.0
nord = 58.0
# Opsætning af rasterfil
kwargs = dict(
driver="GTiff",
height=højde,
width=bredde,
count=antal_bånd,
crs="EPSG:4258", # ETRS89, geografiske koordinater
dtype=np.float64,
transform=rasterio.transform.from_bounds(vest, syd, øst, nord, bredde, højde)
)
with rasterio.open("ny_fil.tif", "w", **kwargs) as grid:
for båndnr in range(1, antal_bånd+1): # bånd er 1-indekserede
data = np.random.rand(højde, bredde)*100
grid.write(data, båndnr)
Læs mere om Rasterio i projektets dokumentation.
Tip
Ved installation af Rasterio følger kommandolinjeværktøjet rio
og med. rio
er et
mere sammenhængende alternativ til GDALs kommandolinjeapplikationer. Læs mere
her.
Fiona¶
Fiona er en moderne pendant til OGR på samme måde som Rasterio er det for GDAL. Fiona tilbyder en mere pythonsk og intuitiv måde at skrive og læse vektordata på.
Bemærk
På engelsk udtales OGR udtales oftest »Ogre« (dansk trold) . Fans af Shrek-serien vil genkende at Fiona er den pæne side af trolden. I Python er Fiona den pæne udgave af OGR.
At læse vektordata med Fiona er lige som at gennemløbe en liste af Python-dicts. Hver feature udgøres af en dict der indeholder information om geometri og attributter. Her et eksempel på en punktfeature fra stednavneregisteret:
{'geometry': {'coordinates': (595507.300561634, 6402053.2234713, -999.0),
'type': 'Point'},
'id': '9480',
'properties': dict([('FEAT_ID', 261996),
('FEAT_KODE', 5010),
('FEAT_TYPE', 'strandpost'),
('AREAL', None),
('SNSOR_ID', 1006651),
('NAVN', 'E320'),
('STAMNAVN', None),
('STAMOPL', 1),
('GODKENDT', '1'),
('INDB_ANTAL', None),
('LAENGDE', None),
('DQ_INDEX', 0),
('TIMEOF_CRE', '2015-07-01'),
('TIMEOF_PUB', '2015-07-01'),
('TIMEOF_REV', None),
('TIMEOF_EXP', '2018-03-31')]),
'type': 'Feature'}
Herunder et eksempel på læsning af punkter fra en shapefil. I eksemplet finder vi det nordligste punkt i datasættet.
import fiona
with fiona.open("eksempel.shp", "r") as punkter:
northing_maks = -1
nordligste_punkt = None
for punkt in punkter:
easting, northing, _ = punkt["geometry"]["coordinates"]
if northing > northing_maks:
northing_maks = northing
nordligste_punkt = punkt
Skrivning af data til en ny fil er selvsagt mere omfattende end at læse data fra en eksisterende fil. I eksemplet herunder oprettes en ny shapefil med afsæt i et schema, der definerer geometritype og attributter for datasættet. Efterfølgende åbnes en ny fil og denne populeres med data.
import fiona
import fiona.crs
SCHEMA = {
"geometry": "Point",
"properties": dict(
[
("navn", "str"),
("kote", "float"),
]
),
}
def ny_feature(navn, kote, easting, northing):
""" Returner ny feature der passer til SCHEMA"""
return {
"geometry": {"type": "Point", "coordinates": (easting, northing)},
"properties": dict(
[
("navn", navn),
("kote", kote),
]
),
}
gnss_stationer = [
ny_feature("SKEJ", 72.028, 573225.054, 6227584.642),
ny_feature("BUDP", 57.890, 719707.285, 6182582.500),
]
with fiona.open(
"ny_fil.shp",
"w",
driver="ESRI Shapefile",
crs=fiona.crs.from_epsg(25832), # ETRS89 / UTM 32N
schema=SCHEMA,
) as stationer:
for station in gnss_stationer:
stationer.write(station)
Tilføjelse af data til et eksisterende datasæt er simplere. Det samme er skrivning af en ny fil med samme opbygning som en eksisterende. Eksempler på disse tilfælde findes i dokumentationen
Tip
Ved installation af Fiona følger kommandolinjeværktøjet fio
og med. fio
er et
mere sammenhængende alternativ til OGRs kommandolinjeapplikationer. Læs mere
her.