Yleistä

Avoimen datan idea, noin yleisellä tasolla, lienee, että yhteisillä rahoilla (~verorahoilla) hankittu tieto pitäisi olla myös yhteisessä käytössä ilman mitään (merkittäviä) rajoituksia tai lisäkustannuksia. Nykyisessä maailmantilanteessa pitää lisäksi tyytyväisenä huomioida, että aivan kaikki verorahat eivät ole (vielä) menneet Ukrainaan tai USA:sta tehtäviin asehankintoihin. Jätetään politikointi kuitenkin (tällä kertaa) tähän ja siirrytään asiaan.

Ilmatieteen laitos

Ilmatieteen laitoksella on tarjolla paljon avointa dataa, tosin sen (laajempi) hyödyntäminen vaatii hieman vaivannäköä. Ei se kuitenkaan mitenkään ylivoimaista ole. Ilmatieteen laitoksen verkkosivuilla on hyvät ohjeet datan lataamiseen ja Githubista löytyy Python rajapinta tiedon hakemiseksi, hyödynnetään sitä.

Auringon säteilyteho

Näin pakkaskelillä ja sähkön hinnan ollessa taivaissa, energia-asiat ovat pinnalla. Katsotaan kevättä odotellessa auringon säteilytehoa (tai sen puutetta) Ilmatieteen laitoksen avoimesta datasta.

Python-ohjelma

Ohjelma toimii Googlen Colab'issa, joka on niitä Googlen vähemmän tunnettuja palveluja. Colab sisältyy kuitenkin samaan (perus)pakettiin, kuin Gmail ja Drive eli se on monilla jo valmiina odottamassa käyttöä. Colabia käytettäessä omalle koneelle ei tarvitse tehdä mitään asennuksia tai muutoksia eli käyttö pitäisi onnistua myös "tietoturvallisella" työkoneella.
Colab on Googlen palvelimessa toimiva Jupyter Notebook, jossa on mm. Python 3-kehitysympäristö.

Python-ohjelma hakee Jyväskylän lentoaseman mittauspisteen tuloksista Auringon kokonaissäteilytehon viimeisen vuorokauden ajalta ja viiden minuutin mittausvälein. Tuloksista piirretään kuvaaja, jonka vaaka-akselilla on aika minuutteina tietojen lataamishetkestä ja pystyakselilla kokonaissäteilyteho.

# !pip install fmiopendata tarvitaan vain ensimmäisellä kerralla
list = !pip list
if len(list.grep("fmiopendata")) == 0:
  !pip install fmiopendata
from IPython.display import clear_output
clear_output()

import datetime as dt
from fmiopendata.wfs import download_stored_query
import re

# Retrieve the latest 24 hour of data
end_time = dt.datetime.utcnow()
start_time = end_time - dt.timedelta(hours=24)
# Convert times to properly formatted strings
start_time = start_time.isoformat(timespec="seconds") + "Z"
end_time = end_time.isoformat(timespec="seconds") + "Z"

sunrad = download_stored_query("fmi::observations::radiation::multipointcoverage", \
                              args=["parameters=GLOB_1MIN", "fmisid=101339", \
                              "starttime=" + start_time, "endtime=" + end_time, \
                              "timestep=5"])
rawdata = sunrad.data

place = "Jyväskylä lentoasema"
meas = "Global radiation"

radiation = []
cntr = 0
sum = 0
min = 9999
max = -9999
for time in rawdata:
  value = rawdata[time][place][meas]['value']
  units = rawdata[time][place][meas]['units']
  radiation.append([time.isoformat(), value, units])

  if value == value: # NaN check
    cntr += 1
    sum += value
    if min > value:
     min = value
    if max < value:
      max = value

radiation.sort()

print("\n{:<20}{:>7}".format("Arvoja (ei Nan): ", cntr))
print("{:<20}{:>7.1f}{}".format("Keskiarvo: ", sum/cntr, " W/m2"))
print("{:<20}{:>7.1f}{}".format("Min: ", min, " W/m2"))
print("{:<20}{:>7.1f}{}\n".format("Max: ", max, " W/m2"))

print("{:<20}{}".format("Alkuaika: ", start_time))
print("{:<20}{}\n".format("Loppuaika: ", end_time))

# Säteilytehon aikajakauma
import matplotlib.pyplot as plt
%matplotlib inline

endtm = dt.datetime.fromisoformat(end_time[0:-1])
time = []
power = []
for meas in radiation:
  tmd = dt.datetime.fromisoformat(meas[0]) - endtm
  tmin = round((tmd.total_seconds())/60, 0) # Aikaero minuutteina
  time.append(tmin)
  power.append(meas[1])

plt.plot(time, power)
plt.show()

Tulokset

Ja tälläiset tulokset saatiin talvisena pakkaspäivänä, kun sitä (lämmitys)energiaa oikeasti tarvittaisiin. Ei tuo nyt kovin lupaavalta näytä. Ainakaan, jos vertaa sitä poliitikkojen puheisiin Suomen tulevaisuudesta (halvan) uusiutuvan energian suurvaltana.

Maanmittauslaitos

Maanmittauslaitos tarjoaa myös avointa dataa, mutta enimmäkseen ladattavina tiedostoina. Kartta-aineistossa reaaliaikavaatimukset lienevät hieman vähäisemmät kuin esim. säätiedoissa ja ladattavat tiedostot ajavat asiansa ihan hyvin.
Kokeillaan korkeustietoa 2m rasterilla ja ASCII Grid-muodossa, joka on selväkielisenä helposti (?) hyödynnettävissä. Korkeusmallin karttalehtijako on 6 x 6 km, joten mittauspisteitä yhdellä karttalehdellä on 9 miljoonaa ja tiedostokoko noin 65 MB. Kun kohdealue sattui vielä melkein karttalehtijaon risteykseen, niin tiedostoja tarvitaan neljä, yhteensä 260 MB! Entä sitten? Sitä vartenhan tietokoneet ovat, että ne voivat käsitellä suuriakin tietomääriä.

Korkeustiedostossa on ensi muutama rivi otsikkotietoja, joissa on esim. vasemman alanurkan paikka, mittauspisteiden jako (2m) ja mittauspisteiden määrä (3000x3000). Otsikkotietojen jälkeen tulee varsinainen korkeustieto alkaen vasemmasta ylänurkasta (länsi-pohjoinen) edeten oikealle kohti itää. Ja sitten seuraava rivi kohti etelää jne. tiedoston loppuun asti.
Paikkatietojen koordinaatisto on suomalainen ETRS TM35FIN, joka on poikittainen Mercatorin-projektio (Gauss-Kruger) keskimeridiaanilla 27ast itäistä pituutta (katsoin tämän Wikipediasta). Tuo kuulostaa jo aika vaikealta, mutta se on vielä helppoa verrattuna siihen laskentaan, jolla WGS84-koordinaatit (esim. GPS ja ~Google Maps) muutetaan ETRS TM35FIN-koordinaateiksi.

Tuo koordinaatistomuunnos vaikuttaa niin vaikealta, että pitää ottaa tekoäly avuksi. (Sainpas mahdutettua "tekoälyn" näillekin sivuille ja nyt sivuston hype-kerroin on kunnossa.) Sen verran perustelua tuohon "tekoäly"-termiin, että koneoppiminen on (pieni) osa tekoälyä ja tarkoitus on käyttää koordinaatistomuunnoksessa opetusdatalla tehtyä koneoppimista.

WGS84 -> ETRS TM35FIN koneoppimisella

Kohdealueena on korkeusmallin karttalehdet M4211G, M4211H, M4213A ja M4213B, joiden keskikohta on 6816000, 332000 (ETRS TM35FIN). Alue on 12 x 12 km, josta otetaan opetusdataksi pisteet tasajaolla kahden kilometrin välein. Yhteensä 49 pistettä (7 x 7). Näille pisteille haetaan vastaavuus WGS84 <-> ETRS TM35FIN esim. Maamittauslaitoksen karttapaikan koordinaatistomuunnoksella.

Koneoppiminen (opetettu):

Lähtötiedoiksi tarvitaan opetusdata.
Arvataan mallin tyyppi:
- Muuttujien ja tulosten määrä. Tässä 2 -> 2.
- Mallin asteluku. Ensimmäisen asteen (lineaarinen) malli olisi yksinkertaisin, mutta pallopinnan kaarevat leveys- ja pituuspiirit aiheuttavat virhettä. Otetaan toisen asteen malli.
- Otetaan laskentakoordinaatiston nollapisteiksi kohdealueen keskikohta. Tämä yksinkertaistaa mallia ja vähentänee (?) virhettä, kun käsitellään vain (suhteellisen) lyhyitä etäisyyksiä ja mallin keskikohta voidaan pakottaa nollaksi (bias/intercept = 0).

Käytetään taas Googlen Colab'ia ja sen Python3-ohjelmointialustaa.
Ohjelman toiminta:

Luetaan opetusdata tiedostosta.
Nollataan opetusdatan keskikohta.
Muodostetaan (koneoppimisella) toisen asteen malli koordinaatistomuunnoksesta.
Tallennetaan laskentamallin vakiot tiedostoon.
Tehdään laskentamallin virhearvio (opetusdatalla).
- Toisen asteen mallilla virhealueeksi saatiin: -0.017...0.015 m, joka on hyvinkin riittävä tähän käyttötarkoitukseen.
- Vertailun vuoksi ensimmäisen asteen mallin virhealue oli n. -10...10 m.


# Koneoppimisella tehty malli koordinaattimuutokseen
# WGS84 -> ETRS-TM35FIN

import os
import numpy as np
import re
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import json

file = "Dataset_HV.csv"
if not os.path.isfile(file):
  !tar -xvf /content/drive/MyDrive/zColab_data/Data.tar

# Opetusdata:
# North, east, latitude, longitude
with open(file, encoding="utf-8") as f:
  f.readline() # Otsikkorivejä ei tarvita
  f.readline()
  lines = f.readlines()

data = []
for line in lines:
  tmp = re.split(";", line)
  data.append([
      float(tmp[0]),
      float(tmp[1]),
      float(tmp[2]),
      float(tmp[3])
  ])
dataset = np.array(data)

# Tekee koordinaattimuunnosmallin ja tallentaa sen json-tiedostoon.
# Keskimmäiset arvot
north_m = np.median(dataset[:,0])
east_m = np.median(dataset[:,1])
for row in dataset:
  if row[0] == north_m:
    if row[1] == east_m:
      ZeroCoords = row
      break

# Nollaus keskikohtaan
dataZero = dataset - ZeroCoords

# Datat mallin opetusta varten
x = dataZero[:,2:4] # Latitude, longitude
y = dataZero[:,0:2] # North, east

# Toisen asteen malli, kaksi muuttujaa, nollapisteen kautta
polyfeatures = PolynomialFeatures(degree=2, include_bias=False)
x_poly = polyfeatures.fit_transform(x)

model = LinearRegression(fit_intercept=False)
model.fit(x_poly, y)

# Tallennetaan mallin laskentavakiot tiedostoon
wgs84_tm35fin = {
    "ZeroCoords": np.ndarray.tolist(ZeroCoords),
    "Coef": np.ndarray.tolist(model.coef_)
    }

with open("wgs84_tm35fin.json", "w", encoding="utf-8") as f:
  f.write(json.dumps(wgs84_tm35fin))

# Virhearvio
test = model.predict(x_poly)
virhe = test - dataZero[:,0:2]
print("\nMallin virhe opetusdatalla:\n")
print("{:3.3f}{}{:3.3f}{}".format(virhe.min(), "...", virhe.max(), " m\n"))

Reitin korkeustieto

Kokeillaan edellä tehtyä laskentamallia ja MML:n korkeusmallitiedostoja reitin korkeusprofiilin muodostamiseen. Lähtötiedoksi piirustellaan Google Maps'issä reitti ja tallennetaan se csv-tiedostoksi. Tiedostossa koordinaatit ovat ~WGS84.
Ohjelman toiminta:

Luetaan reitti- ja korkeusdata-tiedostot.
Muutetaan reitin koordinaatit WGS84 -> ETRS TM35FIN.
Lasketaan reitistä kuljettu matka.
Haetaan reittipisteille korkeustieto.
Piirretään tuloksista korkeusprofiili.

Ohjelman osat:

Lähtötietojen haku
Pääohjelma
- Koordinaatistomuunnos
- Korkeustiedon haku
  - Korkeusdata-tiedoston haku


# Käytetty reittidata on tallennettu Google Mapsistä .csv-muodossa.

import os
import re
import numpy as np
import json
import glob

# Korkeusdata
if not os.path.isfile("M4211G.asc"):
  !tar -xvf /content/drive/MyDrive/zColab_data/Elevations.tar

base_dir = '/content/'
files = glob.glob(base_dir + 'M42*') # Datatiedostot

# Reitti
file = "juoksulenkki.csv"
if not os.path.isfile(file):
  !tar -xvf /content/drive/MyDrive/zColab_data/lenkki.tar

with open(file, encoding="utf-8") as f:
  rawdata = f.readline() # Ensimmäistä riviä ei tarvita
  rawdata = f.readline() # Reitti on toisella rivillä

res = re.search("LINESTRING \((.*)\)\"", rawdata) # LINESTRING sisältö
lines = re.split(",", res.group(1)) # Jaetaan pilkun kohdalta

data = []
for line in lines:
  tmp = re.search("^\s*(.*)", line) # Välilyönnit pois alusta
  tmp = re.split("\s", tmp.group(1)) # Jaetaan välilyönnin kohdalta
  data.append([
      float(tmp[1]),
      float(tmp[0])
  ])

gpsroute = np.array(data) # Reitti: latitude, longitude


# Koordinaatistomuunnos WGS84 -> ETRS TM35FIN
def wgs84_tm35fin(route):
  # Muunnoksen vakioparametrit
  jsonfile = "wgs84_tm35fin.json"
  with open(jsonfile, encoding="utf-8") as f:
    parms = json.loads(f.read())

  Zero = np.array(parms["ZeroCoords"])
  Coef = np.array(parms["Coef"])

  # Laskennan nollapiste
  inp = route - Zero[2:4]

  # Laskennalliset lisäparametrit (toisen asteen malli)
  la = inp[:,0]
  lo = inp[:,1]
  tmp = la ** 2 # Latitude ^2
  inp = np.hstack((inp, np.atleast_2d(tmp).T))
  tmp = la * lo # Latitude * longitude
  inp = np.hstack((inp, np.atleast_2d(tmp).T))
  tmp = lo ** 2 # Longitude ^2
  inp = np.hstack((inp, np.atleast_2d(tmp).T))

  outp = np.matmul(inp, Coef.T)

  # Nollapisteen palautus
  outp = outp + Zero[0:2]

  return outp


# Hakee paikalle datatiedoston
def getFile(point):
  north = point[0]
  east = point[1]
  headers = {}
  for file in files:
    with open(file, encoding="utf-8") as f:
      i = 0
      head = {}
      while i < 6: # Luetaan otsikkotiedot
        line = re.split("\s+", f.readline())
        head[line[0]] = float(line[1])
        i += 1
    headers[file] = head

  datafile = ""
  for file in headers:
    xl = headers[file]["xllcorner"]
    yl = headers[file]["yllcorner"]

    if ((xl <= east) and \
      (east < (xl + headers[file]["ncols"] * headers[file]["cellsize"])) and \
      (yl <= north) and \
      (north < (yl + headers[file]["nrows"] * headers[file]["cellsize"]))):

      datafile = file
      break

  return datafile


# Lukee datatiedoston tiedot
def readData(datafile):
  with open(datafile, encoding="utf-8") as f:
    i = 0
    header = {}
    while i < 6: # Luetaan otsikkotiedot
      line = re.split("\s+", f.readline())
      header[line[0]] = float(line[1])
      i += 1

    ascdata = f.readlines() # Luetaan datarivit

  data = []
  for line in ascdata:
    ascrow = re.split("\s+", line)
    row = ascrow[1:-1] # Alussa ja lopussa on tyhjä
    data.append(row)

  return {"Header": header, "Data": data}


# Hakee paikan korkeuden tiedostosta
def getElevation(route):
  ele = []
  prevdatafile = " "
  for point in route:
    north = point[0]
    east = point[1]
    datafile = getFile(point)
    if datafile != prevdatafile:
      print("Datafile: ", datafile)
      prevdatafile = datafile
      if len(datafile) > 0:
        datadict = readData(datafile)
        header = datadict["Header"]
        data = datadict["Data"]
      else: # Datatiedostoa ei löytynyt
        ele.append(-9999.0)
        continue

    r = (north - header["yllcorner"]) / header["cellsize"]
    r = int(header["nrows"] - r - 1)
    c = int((east - header["xllcorner"]) / header["cellsize"])
    ele.append(float(data[r][c]))

  return ele


# Pääohjelma
# Koordinaatistomuunnos
route = wgs84_tm35fin(gpsroute)

# Matka
delta = route[1:,:] - route[:-1,:]
tmp = delta ** 2
tmp = tmp[:,0] + tmp[:,1]
dist = tmp ** (1/2) # Mittauspisteiden välimatka

odo = [0]
sum = 0
for d in dist:
  sum += d
  odo.append(sum) # Kokonaismatka

# Korkeudet
ele = getElevation(route)

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(odo, ele)
plt.xlabel("Matka / m")
plt.ylabel("Korkeus / m")
plt.show()

Ja tässä tulos lenkkimaastosta.