Van tiff + WGS84 map naar Geotiff in RD formaat

[pxp]

Hoi Ozi gebruikers,

 

ik heb een aantal tiff met OziExplorer map bestanden in WGS84. Deze wil ik graag omzetten naar geotiffs in het Rijksdriehoekstelsel / RD.

 

Ik kan de tiff + tfw (1) bestanden naar geotiff's met WGS84 omzetten via volgende procedure:

 

1) Plaats een bestand met de naam "tfw.dat" in de OziExplorer applicatie directory

2) In OziExplorer --> Options is nu een extra functie om een tfw bestand weg te schrijven

3) FWTools (2) / GDAL gebruiken om met geotifcp.exe een geotiff te maken:

 

geotifcp -e invoer.tfw invoer.tif uitvoer.tif

 

Het bestand uitvoer.tif is nu een geotiff n.a.v. de tif + tfw in WGS84.

 

 

Hoe kan ik deze geotiff in WGS84 eenvoudig omzetten naar het RijksDriehoekStelsel?

 

Of is er een eenvoudigere manier om van een tiff+map in WGS84 formaat een geotiff in RD te maken?

 

Alvast bedankt,

 

Jan

 

 

(1) Het tfw bestand is een informatie bestand met de pixelgrootte en positie om een normale tiff als een geo-gerefereerde tiff in diverse applicaties in te kunnen lezen. Deze zal dan op de correcte coördinaten vallen

(2) FWTools: A library and set of commandline utility applications for reading and writing a variety of geospatial raster (GDAL) and vector (OGR) formats.

1 januari 2006 bewerkt door pxp

[Prof. Y. Lupardi]

Ik heb mij wel eens afgevraagd wat voor soort twf Ozi maakt. Ik heb het nageprobeerd en hetvolgende blijkt:

1.is de kaart geijkt in WGS84 met projectie lat/lon dan is de linkerbovenhoek in decimale graden en de pixelsize ook.

2. is de kaart geijkt in WGS84 met projectie UTM dan gaat ozi over op metrisch en is de linkerbovenhoek in UTM coörds en de pixelsize in meters.

3. is de kaart Rijksdriehoeksmeting geijkt en is de projectie TM met de benodigde parameters om onze RD-projectie te benaderen dan is de twf file in rd-coörds en de pixelsize in meters.

 

Aangezien jij WGS84 geijkt bent zal er conversie op de positie in de twf file noodzakelijk zijn. Met CC40 van RDNAP kan je batchverwerken. Dus als je de twf files uiteenrafelt tot een tabel met posities dan kan je converteren.

Ik zie dat CC40 van UTM ETRF89 (==WGS84) maar ook van decimale graden kan gaan naar RD meters.

 

Maar waarom geotif's in RD?

[pxp]

Bedankt voor de reactie. Of het om geval 1) of 2) gaat qua projectie weet ik niet geheel zeker voor elke map-file.

 

De eerste regels van een voorbeeld bestand zijn :

 

OziExplorer Map Data File Version 2.2

T01.tif

C:\ozi\T01.tif

1 ,Map Code,

WGS 84,WGS 84, 0.0000, 0.0000,WGS 84

Reserved 1

Reserved 2

Magnetic Variation,,,E

Map Projection,Latitude/Longitude,PolyCal,Yes,AutoCalOnly,No,BSBUseWPX,No

 

Ik neem aan dat het in dit bovenstaande geval om een lat/lon projectie gaat.

 

De uitgevoerde tfw vanuit Ozi is:

 

0.00003018

0.00

0.00

-0.00001797

5.45

53.52

 

De laatste 2 regels zie ik bij MMPLL,1, in het map bestand terugkomen, maar de 1e 4 regels zijn waarschijnlijk rekentechnisch bepaald.

 

 

Waarom geotiff,...ehm,...

 

M'n applicatie is gebasseerd op ECW bestanden. Als ik een geotiff heb kan ik deze met GDAL omzetten naar ECW en die vervolgens gebruiken voor m'n eigen navigatiesysteem.

 

Ook is het direct inladen van een geotiff een optie, maar dat gedeelte is nog niet 100% zoals ik het wilde hebben en de ECW bestanden laden zeer snel ook als er meerdere tegelijk aanwezig zijn. Ik ken ECW nog niet zo lang, maar qua snelheid en compactheid van de bestanden is het ideaal.

 

De logica van RD is dat ik de ontvangen NMEA posities realtime converteer naar RD. Dit is voor m'n doeleinden idealer. Ik heb de beschikking over een aantal esri shape bestanden die in RD zijn en ben met die logica gestart. Ook heb ik meer gevoel bij een RD coordinaat dan bij een WGS84 coordinaat. Zal wel iets persoonlijks zijn

 

Ik heb een aantal OziExplorer jpeg en tiff bestanden met de Ozi map-files in WGS84.

 

Als ik de map-files batchgewijs kan omzetten naar RD zou het ook al een stuk helpen, althans als de tfw bestanden dan correct zijn om een RD geotiff aan te maken.

 

'k hoop dat deze uitleg misschien een alternatief bied boven het letterlijk converteren van elke map-bestand.

[Prof. Y. Lupardi]

Inderdaad zijn jouw tfw files niet metrisch. De linkerbovenhoek positie is in decimale graden.

Maar ook de celsize (pixelgrootte) is in dec. graden. Ik zie dat je pixels niet vierkant zijn maar ik bedenk net dat graden in de lengte en in de breedte anders uitmeten.

De pixelsize berekent Ozi met een 'save mapfile' steeds opnieuw en die staat bij

MM1B, [decimale meters]

De kaarthoeken (die je met de hand kan plaatsen als de kaart een witte rand heeft) staan in MMPLL,1, [dec. graden WGS84], [dec. graden WGS84] en worden ook bij een save herberekend.

Het monnikenwerk zou neer komen op:

1. lezen van da map-file

2. linkerbovenhoekpositie (0,0) in CC40 typen, converteren.

3. pixelgroote lezen

3. een twf-file maken en de pixelgroote plus de RD coörds intypen.

etc.

 

Ik zie het niet zo gauw maar kan Ozi slechts een tfw maken voor de huidig ingeladen kaart? Dan is tfw maken tijdrovend.

Wat valt er te automatiseren (=programmeren, macro schrijven e.d.) ?

Een tfw file is wat simpeler uiteen te pluizen dan een mapfile. Minder programmeerwerk.

CC40 kan je een tabel met posities voederen ter conversie.

Vanuit de nieuwe tabel (RD) weer tfw files maken.

 

Als je grootschalige kaarten hebt dan is de projectie niet zo vreselijk belangrijk maar de mapdatum wel. Zou je de mapfiles niet lat/lon geven maar UTM dan zal Ozi metrische UTM (WGS84) getallen in de tfw file zetten. Daar kunnen de meeste GIS pakketten wel tegen.

 

Je enthousiasme voor ECW kan ik volledig delen. Zeker nu Ozi CE ook ze kan gebruiken (mits je CPU van het ARM-type is) .

Er is een progje ECW header editor waarin ik de geoinfo in de ECW file kan zetten. (Met de hand).

Met ERMapper kan ik batchgewijs ECW(of gif, jpg, tif) + tfw omzetten tot ECW (geo).

 

Zitten de honderden map files een beetje regelmatig in elkaar dan zou ik wat gaan programmeren ervoor in m-code (MatLab) , even mijn progjes voor de berekeneningen WGS84 -> RD opzoeken en vervolgens ers-files gaan genereren voor ERMapper.

Dan de ers-files batchgewijs door ERMapper halen die dan daarvan samen met de tiff's de geo-ecw's maakt.

Maar nu wordt het wel heel erg technisch. Misschien per PM de verdere detailprobleempjes uitwerken?

[pxp]

Bedankt voor de erg uitvoerige reactie, perfect

 

Ik was al bang dat het geen eenvoudige klus was

 

ERmapper heb ik zelf nog niet voor het converteren gebruikt, dus 'k zal me eens verdiepen in ers-files.

 

Ik zal je een PM sturen. Doe ik dat ook een keer

[Gast Keessie]

Via google hier terecht gekomen. Ook ik ben bezig om Oziexplorer map files te converteren naar geotiff, maar weet niet of je verder ben gekomen.

 

In mijn geval heb ik heb zowel van beide situaties werkende voorbeelden, dus ik heb ozi-mapfiles + ogf2 en worldfile + tif. Om beide situaties in 1 pakket te laten draaien moet ik de oziexplorer map files converteren naar twf bestanden. Beide syntax-en zijn niet heel erg lastig te begrijpen het enige is de conversie naar de coordinaten stelsel.

 

Ik weet niet hoe goed je programmeer kunsten zijn, maar in mijn situatie komt het neer op het converteren van ogf2 naar TIF (wellicht eerst via BMP) en vervolgens de WGS 84 coordinaten uit de map-file halen en converteren naar RD- twffile.

 

Uit eerdere projecten heb ik nog de volgende conversie code, weliswaar in PHP maar dat moet geen probleem vormen.

 

function rd_nb_ol_x_y($nb,$ol)
{
  global $x_r,$y_r,$gods_mode;

$f=$nb;
$l=$ol;

$x0  = 155000.00;
$y0  = 463000.00;

$f0 = 52.15616056;
$l0 =  5.38763889;

$c01=190066.98903;  $d10=309020.31810;
$c11=-11830.85831;  $d02=  3638.36193;
$c21=  -114.19754;  $d12=  -157.95222;
$c03=   -32.38360;  $d20=    72.97141;
$c31=    -2.34078;  $d30=    59.79734;
$c13=    -0.60639;  $d22=    -6.43481;
$c23=     0.15774;  $d04=     0.09351;
$c41=    -0.04158;  $d32=    -0.07379;
$c05=    -0.00661; $d14=    -0.05419;
$d40=    -0.03444;
$df=($f -$f0) * 0.36;
$dl=($l - $l0) * 0.36;


$dx =$c01*$dl + $c11*$df*$dl + $c21*pow($df,2)*$dl + $c03*pow($dl,3);
$dx=$dx+$c31*pow($df,3)*$dl + $c13*$df*pow($dl,3) + $c23*pow($df,2)*pow($dl,3);
$dx=$dx+$c41*pow($df,4)*$dl + $c05*pow($dl,5);
$x=$x0 + $dx;
$x_r=round(100*$x)/100;


$dy =$d10*$df + $d20*pow($df,2) + $d02*pow($dl,2) + $d12*$df*pow($dl,2);
$dy=$dy+$d30*pow($df,3) + $d22*pow($df,2)*pow($dl,2) + $d40*pow($df,4);
$dy=$dy+$d04*pow($dl,4) + $d32*pow($df,3)*pow($dl,2) + $d14*$df*pow($dl,4);
$y=$y0 + $dy;
$y_r=round(100*$y)/100;

if($debug=='1'){echo"<br>x wordt :x=$x_r<br>y wordt y=$y_r<br>";}
return;

}

 

Laat me eens weten of je (jullie) verder bent gekomen..

[pxp]

Hoi Keessie,

 

bedankt voor je info. In principe is het probleem inderdaad om de Ozi map bestanden om te zetten naar tfw.

 

Probleem is dat de projectie WGS84 is met meerdere ijkpunten. Ik heb wel een routine om WGS84 om te zetten naar RD en andersom, maar probleem is dus de meervoudige ijkpunten die in de map-file gedefineerd zijn op willekeurige posities. Het lijkt me een affine transformatie waar Ozi haar bestand mee refereerd en dat is het lastige in dit geval.

 

Dus is het denk ik niet een questie van alleen maar de coordinaten vertalen van WGS84 naar RD, maar een 'image transformatie'. Dus geen translatie [verschuiving] maar een transformatie [verplaatsen/roteren/uitrekken]....toch?

 

Met welke soorten bestanden is het je gelukt om de conversie toe te passen?

 

Hoor graag van je,

 

Jan

[Prof. Y. Lupardi]

Jouw probleem is dat de ijkpunten overal kunnen liggen en je moet de kaarthoeken in RD hebben.

Zoals al gezegd: met een 'save mapfile' berekent Ozi zaken opnieuw.

Onderaan de regels met PMPLL, bevatten de decimale WGS84 graden van de geplaatste kaarthoeken. De positie van de kaarthoeken staat in pixels in MMPXY, en de rechter onderhoek van de kaart heb je te pakken met de onderste regel IWH, Map Image....

Nu slechts de kaarthoekposities omzetten naar RD.

 

"Probleem is dat de projectie WGS84 is met meerdere ijkpunten. " moet zijn:

"Probleem is dat de mapdatum WGS84 is met meerdere ijkpunten. "

Verwar projectie niet met mapdatum.

Een reprojection van een kaart is alleen nodig als die b.v. als ondergrond gebruikt wordt voor lagen met andere gegevens die in een andere projectie staan. Of als je kaartdelen aaneen moet voegen die verschillende projecties hebben.

 

Jouw kaarten hebben lat/lon als projectie oftewel: onbekende projectie!

[pxp]

Hoi Prof,

 

klopt. Blijft toch lastig praten, projectie etc... maar ik bedoel het goed

 

Heb het een tijdje laten liggen,... dit weekend eens kijken wat er mogelijk is.

[pxp]

Ben er nog weer eens mee aan het stoeien, maar ik denk dat de tiff via de map-file door Ozi intern geprojecteerd word middels de calibratiepunten.

 

Wanneer ik een map-file save dan heb ik o.a.:

 

MM0,Yes

MMPNUM,4

MMPXY,1,0,0

MMPXY,2,26229,0

MMPXY,3,26229,20633

MMPXY,4,0,20633

MMPLL,1, 1.086037, 53.771922

MMPLL,2, 8.252652, 53.790787

MMPLL,3, 8.027491, 50.449450

MMPLL,4, 1.356293, 50.432157

MM1B,18.036635

LL Grid Setup

LLGRID,No,No Grid,Yes,255,16711680,0,Geen label,0,16777215,7,1,Yes,x

Other Grid Setup

GRGRID,No,No Grid,Yes,255,16711680,Geen label,0,16777215,8,1,Yes,No,No,x

MOP,Map Open Position,0,0

IWH,Map Image Width/Height,26229,20633

 

Als ik het goed begrijp is:

 

* linkerboven pixelpositie: MMPXY,1,0,0 [0,0 pixelcoord]

* linkerboven l/l positie: MMPLL,1, 1.086037, 53.771922

* pixelgrootte: MM1B,18.036635

 

Voor het maken van de tfw gebruik ik dus de linkerbovenpositie [1.086037, 53.771922] en de pixelgrootte

 

In RD is deze positie: -128499.89 , 651416.46

 

 

In Ozi is de positie van de kaart in combinatie met een GPS coord goed. Wanneer ik van de map file's gegevens een tfw maak met linksboven-punt-in-RD-coordinaten en de pixelgrootte en dit vervolgens test in een RD-applicatie dan is de positie van de raster foutief. Zit er zo'n 600 meter naast.

 

De tiff [origineel] zal uitgerekt moeten worden denk ik en kan niet zomaar met de twf naar RD worden omgezet.

[Prof. Y. Lupardi]

Het betreft hier een kleinschalige kaart. Die loopt over een groot bereik van 1 tot 8 west en 50 tot 54 Noord. Ver buiten het bereik van onze RD.

Ook UTM gaat hier niet omdat je hier over diverse zones heel loopt.

M.a.w. er is geen metrisch raster mogelijk voor deze kaart.

Wat inhoud dat een tfw file onbruikbaar is. Die moet metrische (of desnoods surveyers feet) eenheden hebben.

Een tfw file die posities en afstanden in graden bevat gaat niet. Positie linker bovenhoek, pixelsize in x-richting, pixel size in y-richting.

Maar bij graden is de horizontale afstand aan de bovenzijde van de kaart kleiner dan aan de onderzijde. (of lopen je N/Z lijnen parellel?)

Dit gaat ook op voor geotif: die verwacht coördinaten in een metrisch rechthoekig raster.

Je schreef:

"M'n applicatie is gebasseerd op ECW bestanden. Als ik een geotiff heb kan ik deze met GDAL omzetten naar ECW en die vervolgens gebruiken voor m'n eigen navigatiesysteem."

Mijn ECW Image Header Editor geeft hetvolgende:

Projection Detail: Datum, Projection, Units (meters, deg, unknown),

Registration Details: Top left X coordinate, top left Y-coordinate, Cell witdth X cell height Y

 

Heb je enig idee of de projectie de eigenschap 'afstanden' bewaard heeft?

[pxp]

even een snelle tip tussendoor,...

 

ERDAS ViewFinder 2.1

 

"one of the few (if not the only) free GIS image viewer that can also enhance and reproject raster data."

 

 

download links:

ftp://ftp.gi.leica-geosystems.com/softwar...VF_Setup_21.EXE

http://www.ermapper.com/downloads/download..._VERSION_ID=151

 

Er is dus een extra ecw plugin beschikbaar [2e link].

 

Her-projectie Voorbeeld:

http://www.digitalgrove.net/image_reprojec..._with_erdas.htm

 

Mogelijke projecties staan in deze lijst: http://www.digitalgrove.net/images/erdas_proj.gif

25 februari 2006 bewerkt door pxp

[pxp]

Ik heb het ook met de header editor tool ingetikt en als units gebruik ik meters.

Dit is ook zo'n 1000 meter fout

25 februari 2006 bewerkt door pxp

[Prof. Y. Lupardi]

Dank voor de aanwijzing betreffende de plugin voor ERDAS Viewfinder. Die stond bij mij te boek als voor MrSid spullen only"." Nu niet meer. Ik kan nu ook .ers files inladen die doorlinken naar ecw plaatjes.

 

PS. Als je als units 'degrees' geeft en alles in decimale graden opgeeft?

[Gast Keessie]

Ik ben ergens afgehaakt, hou het draadje nog in de gaten en zal mn voortgang qua conversie nog posten.

[pxp]

Toevallig nog [free] software gevonden om naar RD te transformeren?

Of naar ETRS

 

ViewFinder heeft namelijk geen RD parameters.