Projectie van KNMI radarbeelden
Iwan Holleman / Hans Beekhuis
Inhoud:
- Geografische positie van de radars
- Geografische positie naar radarbeeld
- Radarbeeld naar Geografische positie
- Projectie parameters van het radarbeeld met 1.0km x 1.0km pixels
- Projectie parameters van het radarbeeld met 2.5km x 2.5km pixels
- Projectie berekeningen mbv Proj.4
Geografische positie van de radars:
De posities van de twee weerradars van het KNMI zijn:
| Radar | Lengtegraad | Breedtegraad | Lengtegraad | Breedtegraad |
| De Bilt | 5o10'42.04"E | 52o06'06.04"N | 5.17834E | 52.10168N |
| Den Helder | 4o47'23.90"E | 52o57'12.02"N | 4.78997E | 52.95334N |
Geografische positie naar radarbeeld:
De KNMI radarbeelden hebben een polair stereografische projectie. Een stereografische projectie gaat uit van een centraal punt (oorsprong) en een referentie meridiaan, in ons geval de noordpool en de Greenwich meridiaan. Een geografische positie met lengtegraad (L) en breedtegraad (B) wordt omgerekend naar een afstand tot de oorsprong en een hoek ten opzichte van de referentie meridiaan (L0). Voor een stereografische projectie geldt dat hoeken behouden zijn (conforme projectie). Meridianen worden rechte lijnen vanuit de noordpool en breedtecirkels worden cirkel gecentreerd rond de noordpool. De projectie wordt uitgevoerd op een ellipsoide aarde met een straal op de evenaar van Re en een straal op de polen van Rp. Voor de omrekening van de breedtegraad naar de projectie afstand R(B) geldt:
Eccentriciteit e wordt berekend uit de aardstralen (e=0 --> sferisch):
- e = sqrt( 1 - (Rp/Re)^2 )
De ellipsoide factor F(B) volgt uit:
- f(B) = (1-e)*(1+e*sin(B)) / [ (1+e)*(1-e*sin(B)) ]
- F(B) = f(B)^(e/2) / sqrt(1-e^2)
De projectie afstand R(B) wordt nu:
- R(B) = 2 * Re * tan(45-B/2) * F(B)
De projectie afstand wordt geschaald zodanig dat de schaalfactor gelijk is aan 1 op de breedtegraad van ware schaling (Bw). Voor de schaalfactor Z(Bw) geldt:
De ellipsoide schaal factor G(Bw) volgt uit:
- g(Bw) = (1-e)*(1+e*sin(Bw)) / [ (1+e)*(1-e*sin(Bw)) ]
- G(Bw) = g(Bw)^(e/2) * sqrt[ (1-(e*sin(Bw))^2) / (1-e^2) ]
De schaal factor Z(Bw) is nu gelijk aan:
- Z(Bw) = 2 * tan(45-Bw/2) * G(Bw) / cos(Bw)
De Cartesische afstanden worden berekend in een stelsel waarbij de y-as parallel ligt aan de referentie meridiaan en naar het noorden wijst en de x-as loodrecht daarop. De Cartesische afstanden X tot de referentie meridiaan en Y tot de oorsprong van de projectie zijn:
- X(L,B) = R(B) * sin(L-L0) / Z(Bw)
- Y(L,B) = -R(B) * cos(L-L0) / Z(Bw)
Tenslotte worden de Cartesische afstanden omgerekend naar beeldpunt coordinaten I en J met behulp van de beeldpuntgrootte (PX en PY) en de verschuiving van het eerste beeldpunt (OI en OJ):
- I(L,B) = X(L,B) / PX - OI
- J(L,B) = Y(L,B) / PY - OJ
Deze beeldpunt coordinaten kunnen worden omgezet naar beeldpunt nummers door een truncatie (int) operatie.
Radarbeeld naar Geografische positie:
Voor de omrekening van beeldpunt coordinaten of nummers naar geografische positie moet bovenstaande procedure in omgekeerde volgorde worden uitgevoerd. Bij 1 berekeningsstap is een iteratie stap vereist. Beeldpunt nummers moet eerste worden omgerekend naar beeldpunt coordinaten door 0.5 erbij te tellen. Vervolgens kunnen uit de beeldpunt coordinaten de Cartesische projectie afstanden worden berekend:
- X(I,J) = (I + OI) * PX
- Y(I,J) = (J + OJ) * PY
De afstand tot de projectie oorsprong en de hoek met de lengtegraad volgen uit:
- R(I,J) = Z(Bw) * sqrt( X(I,J)^2 + Y(I,J)^2 )
- L(I,J) = L0 + atan2( X(I,J) , -Y(I,J) )
De breedtegraad kan worden berekend uit de projectie afstand door een twee staps iteratie. Tijdens de eerste stap wordt bij de inversie van de projectie afstand het "ellipsoide deel" F(B) verwaarloosd. Met behulp van deze eerste schatting van de breedtegraad B1 kan de ellipsoide correctie factor worden berekend. Vervolgens kan een nieuwe schatting van de breedtegraad B worden berekend:
- B1 = 90 - 2 * atan( R(I,J) / (2*Re) )
- B = 90 - 2* atan( R(I,J) / (2*Re*F(B1)) )
Projectie parameters van 1x1 km radarbeeld:
De parameters van het KNMI radarbeeld staan weergegeven in onderstaande tabel. Met behulp van bovenstaande formules kan voor elk beeldpunt een geografische positie worden berekend of voor elke geografische positie het bijbehorende beeldpunt worden gevonden. Het Aardemodel dat hierbij wordt gehanteerd is WSG84.
| Parameter | Waarde |
| Projectie | Stereografisch |
| Oorsprong van projectie | L0=0.0E en B0=90.0N |
| Breedtegraad met ware schaal | Bw=60.0N |
| Straal op evenaar en pool (WGS-84) | Re=6378.137 km en Rp=6356.752 km |
| Verschuiving van beeld | OI=0 en OJ=3650.0 |
| Grootte van beeldpunten | PX=1.0 km en PY=-1.0 km |
| Aantal beeldrijen en kolommen | 765 en 700 |
Met deze parameters kunnen bijvoorbeeld de hoekpunten van het radarbeeld worden berekend:
| Hoekpunt | Lengtegraad | Breedtegraad | XXX[km] | YYY[km] |
| Noordwest | 0.000E | 55.974N | 0.000 | -3650.000 |
| Noordoost | 10.856E | 55.389N | 700.000 | -3650.000 |
| Zuidoost | 9.009E | 48.895N | 700.000 | -4415.000 |
| Zuidwest | 0.000E | 49.362N | 0.000 | -4415.000 |
Projectie parameters van het 2.5 x 2.5 km radarbeeld:
(N.B. dit product is vervangen door de 1km x 1km producten) Het aarde model dat gebruikt is niet wsg84 maar het model dat door het "rijksdriehoek amersfoort stelsel" werd gebruikt. De parameters van het KNMI radarbeeld staan weergegeven in onderstaande tabel. Met behulp van bovenstaande formules kan voor elk beeldpunt een geografische positie worden berekend of voor elke geografische positie het bijbehorende beeldpunt worden gevonden.
| Parameter | Waarde |
| Projectie | Stereografisch |
| Oorsprong van projectie | L0=0.0E en B0=90.0N |
| Breedtegraad met ware schaal | Bw=60.0N |
| Straal op evenaar en pool | Re=6378.388 km en Rp=6356.912 km |
| Verschuiving van beeld | OI=0 en OJ=1490.906 |
| Grootte van beeldpunten | PX=2.5 km en PY=-2.5 km |
| Aantal beeldrijen en kolommen | 256 en 256 |
Met deze parameters kunnen bijvoorbeeld de hoekpunten van het radarbeeld worden berekend:
| Hoekpunt | Lengtegraad | Breedtegraad | XXX[km] | YYY[km] |
| Noordwest | 0.000E | 55.296N | 0.0 | -3727.291 |
| Noordoost | 9.743E | 54.818N | 640.0 | -3727.291 |
| Zuidoost | 8.337E | 49.373N | 640.0 | -4367.291 |
| Zuidwest | 0.000E | 49.769N | 0.0 | -4367.291 |
Projectie berekeningen mbv Proj.4:
De projectie berekeningen kunnen ook worden uitgevoerd met behulp van de "proj.4" library die door de USGS is ontwikkeld en wereldwijd in heel veel applicaties wordt gebruikt. Meer informatie over de proj.4 library kan worden verkregen via USGS Proj.4 Projection Library
De projectie van de KNMI radarbeelden kan worden uitgevoerd met behulp van de volgende proj.4 string die de complete geografische proejctie beschrijft:
"+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.388 +b=6356.906"
Terug naar de index
Hans Beekhuis