Polynomial transformation / Georeferencing - de wiskunde!?

[tonyisgaaf]

Ik heb een widget, "NL Weerradar", voor Mac OS X Dashboard. Daarin worden een aantal weerkaarten van verschillende bronnen geladen. Nu speel ik al een tijd met het idee een functie toe te voegen waarmee je een "pin" op een kaart kan plaatsen aan de hand van het opgeven van een lokatie. De latitude/longitude wordt dan opgezocht m.b.v. een geocoding service en vertaald naar het pixel grid van de weerkaart image. Voor deze vertaalslag zijn parameters nodig.

Het "reverse referencing" van het pixel grid van een bitmap image naar de corresponderende latitude/longitude heet "georeferencing" en het wiskundige principe heet "polynomial transformation", ben ik na enig zoeken achter. Bijvoorbeeld Mathlab heeft een dergelijke functie ingebouwd.
De werking is simpelweg gebaseerd op het bekend zijn van bepaalde punten op de kaart en de pixel waarden laten corresponderen met de latitude/longitude waarden. Hieruit volgt een stelsel van vergelijkingen, wat uiteraard valt op te lossen bij voldoende referencing points.

De meest simpele vorm van "georeferencing" kan d.m.v. een "linear transformation" met de vorm:

code:

x = a*longitude + b
y = c*latitude + d

maar dit geldt alleen wanneer de bitmap kaart qua Noord-Zuid oriëntatie klopt en de latitude en longitude gelijkvallen met de horizontale en verticale pixel-rijen in de bitmap.

Dus valt de transformatie te benaderen met een tweedegraads polynoom (volgens mij):

code:

x = a*lon^2 + b*lon + c
y = d*lat^2 + e*lat + f

Daarnaast kunnen de kaarten zijn geroteerd, in het algemeen geldt dan:

code:

x' = x*cos(alpha) - y*sin(alpha)
y' = x*sin(alpha) + y*cos(alpha)

Substitutie van x en y in x' en y' geeft dan:

code:

x' = (a*lon^2 + b*lon + c)*cos(alpha) - (d*lat^2 + e*lat + f)*sin(alpha)
y' = (a*lon^2 + b*lon + c)*sin(alpha) + (d*lat^2 + e*lat + f)*cos(alpha)

Dit stelsel met zeven onbekenden (a, b, c, d, e, f, alpha) heb ik als test m.b.v. een tiental bekende punten (x, y, longitude, latitude) in Maple ingevoerd en in verschillende combinaties gesolved, maar daar komen elke keer (in verschillende combinaties) zulke verschillende combinaties uit (totaal verschillende hoeken bijvoorbeeld), dat ik niet goed weet hoe ik dit stelsel kan oplossen en daarnaast met behulp van de niet gebruikte punten de benadering kan verbeteren.

Verder lees ik op de pagina's van een aantal tools die dit kunnen, zoals GCPWorks => Polynomial transformations dat er zeven punten nodig zijn voor het oplossen van een tweedegraads polynomial transformation, maar voor het stelsel vergelijkingen heb je ogenschijnlijk aan vier punten (= acht vergelijkingen: 4 maal x, 4 maal y) genoeg.

Weet iemand of de vergelijkingen kloppen en hoe ik het georeferencen kan verbeteren? Als er behoefte aan is kan ik de Maple worksheet en/of de tabel met pixel grid en lat/lon waarden online zetten.

[Towap]

Hoi

Ik ben geen programmeur maar transformaties is wel een beetje mijn ding. Volgens mij maak je het wat te moeilijk. De bitmap kan afwijken van de lon en lat qua scalering en rotatie. (Tenzij je ervan uit gaat dat er ook een ongelijke schaling van de assen mogelijk is maar ik denk dat dat niet de bedoeling is)

concreet corrigeert de eerste set vergelijkingen voor herschalen(a) en nieuwe oorsprong (b):

code:

x=a*lon + b 
y=a*Iat + c

a is in beide situaties gelijk tenzij de x en y as dus apart kunnen geschaald worden: je gebruikt dus schaal a en nieuwe oorsprong (b,c)

Nu de rotatie, deze kan inderdaad voorkomen en wordt inderdaad beschreven door:

code:

x' = x*cos(alpha) - y*sin(alpha)
y' = x*sin(alpha) + y*cos(alpha)

In totaal zijn er dus vier onbekenden: alpha, a, b en c. Het stelsel is ook lineair en dus oplosbaar. Merk op dat je alpha niet hoeft te berekenen: sin(alpha) of cos(alpha) zijn onderling verbonden via de regel dat de som van de kwadraten van de twee 1 is. Je hoeft dus enkel cos(alpha) of sinus(alpha) uit te rekenen.

Met matrices is het nog veel simpler, maar tenzij je matlab gebruikt denk ik niet dat je daar veel mee bent.

[tonyisgaaf]

Towap schreef op donderdag 16 juli 2009 @ 18:18:
Hoi

Ik ben geen programmeur maar transformaties is wel een beetje mijn ding. Volgens mij maak je het wat te moeilijk. De bitmap kan afwijken van de lon en lat qua scalering en rotatie. (Tenzij je ervan uit gaat dat er ook een ongelijke schaling van de assen mogelijk is maar ik denk dat dat niet de bedoeling is)
[...]

Misschien was ik daarin niet duidelijk, maar dat is juist het probleem: de testafbeelding is (wanneer handmatig teruggeroteerd en geschaald) niet gelijk aan een andere afbeelding waarvan ik "zeker" weet dat deze "conformal" geprojecteerd is (zie voor definitie "conformal"). Daarom gebruik ik dus ook deze ingewikkelde variant, omdat ik zeker weet dat de schaal in beide richtingen anders is.

De lineaire transformatie had ik al geprobeerd (wel met a, b, c, d en alpha) en dan klopt de uitkomst voor de ingevoerde punten P en Q, maar voor controlepunten die afwijken van de lijn P-Q loopt de afwijking snel op.

[pedorus]

Gaat het niet gewoon om een Mercatorprojectie? Aangezien Buienradar ook een Google Maps versie heeft, neem ik aan dat ze die gebruiken...

[tonyisgaaf]

pedorus schreef op donderdag 16 juli 2009 @ 20:01:
Gaat het niet gewoon om een Mercatorprojectie? Aangezien Buienradar ook een Google Maps versie heeft, neem ik aan dat ze die gebruiken...

Het klopt inderdaad dat Google Maps een Mercator projectie is. Als ik echter de standaard Buienradar afbeelding over de Google Maps afbeelding van (ongeveer) dezelfde schaal leg, krijg ik het niet voor elkaar de rotatie en de schaal te laten overeenkomen. Als ik bijvoorbeeld Cambridge en Koblenz (zo'n beetje uitersten op de Buienradar kaart) van beide kaarten over elkaar laat vallen, dan vallen andere plaatsen (zelfs Norwich vlak bij Cambridge) al niet meer samen.
De Buienradar kaart gebruikt dus een andere projectie en is daarnaast geroteerd. Met behulp van de tweedegraads polynoom transformatie had ik gehoopt een benadering daarvan te kunnen krijgen.

Ik snap nog steeds niet waarom andere bronnen zeggen dat er minimaal zeven punten nodig zijn voor het oplossen van het stelsel, terwijl er naar mijn idee vier genoeg zijn (voor acht vergelijkingen). Ik heb sterk het idee dat ik iets over het hoofd zie.

Even een afbeelding als voorbeeld, Buienradar afbeelding over Google Maps (terrain):
(klik voor groter)

[ Voor 9% gewijzigd door tonyisgaaf op 17-07-2009 00:18 ]

[pedorus]

Hm, ja. wap.buienradar.nl maakt gebruik van proj. Stuur anders even een mailtje naar hem welke transformatie(s) hij gebruikt? Ik ben wel benieuwd welke projectie ze dan gebruiken...