Re: [FreeGeoCZ] Konverze z S-JTSK pomoci gdalwarp

[Michal <michal AT lightcomp.cz>]

Ahoj,

tvuj zpusob vypada mnohem elegantneji, ale bohuzel mi nejak nefunguje. 
Zkousel jsem:

gdal_translate -of vrt -a_srs '+proj=krovak +lat_0=49.5 
+lon_0=24.83333333333333 +alpha=0 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel 
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m +no_defs' \
  -a_ullr -602500 -1034000 -593500 -1043000 \
  -expand rgba \
  tmp.tif pokus.vrt

gdal2tiles.py pokus.vrt

Ale vysledek je uplne mimo, BoundingBox v tilemapresource.xml je:
<BoundingBox minx="-9.32794713297028" miny="-5.41234958682012" 
maxx="-9.24815873556103" maxy="-5.33150121124936"/>

Jinak jsem zjistil, ze moje worklow funguje dobre, pokud konvertuju pres 
longlat wgs84. T.j. ne primo S-JTSK -> TMS ale nejdrive S-JTSK -> 
"+proj=longlat +datum=WGS84" a potom "+proj=longlat +datum=WGS84" -> 
TMS. Mimochodem, pokud je to znama vec jak pises, ze S-JTSK->TMS se 
neprovede spravne, nemel by gdalwarp tuto transformaci odmitnout nebo 
alespon zarvat?

GeoTIFF v "+proj=longlat +datum=WGS84" mi gdal2tiles nechce zpracovat a 
musim konvertovat do TMS. To ma nejaky duvod, nebo je to bug?

Take jsem zkousel pouzit pri transformaci grid 
(http://grass.fsv.cvut.cz/gwiki/S-JTSK-Grid), ale zatimco cs2cs funguje 
spravne s "+proj=krovak +ellps=bessel +nadgrids=czech", tak gdalwarp mi 
vysledek hazel zase nekam do more a fungovalo to jen s "+proj=krovak 
+lat_0=49.5 +lon_0=24.83333333333333 +alpha=0 +k=0.9999 +x_0=0 +y_0=0 
+ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m 
+nadgrids=czech +no_defs". Jenom si nejsem jistej, zda se v tomto 
pripade grid pouziva spravne.

Michal

Klokan Petr Přidal wrote:
> Ahoj,
>
> nezkousel jsem tvoje workflow ale obecne bych pouzil radeji postup
> podobny tomuto:
>
> gdal_translate -of vrt -a_srs '+proj=krovak +lat_0=49.5
> +lon_0=24.83333333333333 +alpha=0 +k=0.9999 +x_0=0 +y_0=0
> +ellps=bessel +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
> +units=m +no_defs' \
> -a_ullr -602501 -1034000 -602500 -1034001
> pokus_1034001-1034000_602501-602500.png pokus.vrt
>
> V tento moment by mel pokus.vrt vracet stejne vysledky jako cs2cs) a
> raster je korektne georeferencovan do S-JTSK.
>
> a pak primo:
>
> gdal2tiles.py pokus.vrt
>
> Proc?
> Gdal2tiles bude warpovat do cilove projekce sam tak jako tak... Pokud
> warpujes do spherical mercatoru pomoci proj.4 definice do vysledku se
> ulozi WKT bez informace o nulove +nadgrids transformaci, spravne by
> tam melo byt toto WKT:
> http://www.spatialreference.org/ref/sr-org/6/prettywkt/
> V gdalinfo tedy chybi EXTENSION parameter...
>
> Proto je vysledek a dalsi zpracovani podobneho rastru
> problematictejsi, protoze se mercator provadi na skutecnem WGS84
> elipsoidu namisto na kouli - to zpusobi transformacni posun ktery
> vidis.
>
> Jestli je nejaka dalsi chyba na urovni S-JTSK reference vstupnich dat
> na strane CENIE, jak pise Štěpán Kafka, to nevim... nezkousel jsem.
>
> Jinak GDAL umi cist i WMS takze gdal2tiles by melo byt schopne
> vytvaret dlazdice primo z dotazu na geoportal.cenia.cz (stejne jako z
> jakehokoliv jineho WMS serveru) ale musi se pro takove WMS vytvorit
> XMLko popisujici okno, url serveru, projekci a velikost jednotlivych
> dotazu...
>
> Klokan
>
> 2009/4/2 Michal <michal AT lightcomp.cz>:
> > Štěpán Kafka wrote:
> > > NO, obávám se, že něco posunutého mají na CENII, taky jsme s tím měli
> > > potíže
> > To jsem si myslel puvodne taky, ale lidi s uspechem tyhle mapy pouzivaji v
> > OziExploreru a souradnice rohu v map souboru odpovidaji vystupu z cs2cs.
> > Take jsem udelal dalsi primitivni pokus s obrazkem velikosti 2x2 pixelu a
> > 1mx1m v S-JTSK:
> >
> > convert pokus_1034001-1034000_602501-602500.png tmp.tif
> >
> > gdal_translate \
> >  -a_srs '+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=0
> > +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
> > +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m +no_defs' \
> >  -gcp 0 0 -602501 -1034000 \
> >  -gcp 0 1 -602501 -1034001 \
> >  -gcp 1 0 -602500 -1034000 \
> >  -gcp 1 1 -602500 -1034001 \
> >  tmp.tif tmp_ref.tif
> >
> > gdalwarp \
> >  -s_srs '+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=0
> > +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
> > +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m +no_defs' \
> >  -t_srs '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
> > +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs' \
> >  tmp_ref.tif pokus_1034001-1034000_602501-602500.tif
> >
> > podle vypoctu cs2sc jsou rohy nasledujici:
> > $ cs2cs -f "%.10f" +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333
> > +alpha=0 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel
> > +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +units=m +no_defs +to
> > +proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs
> > -602501 -1034000
> > 16.3561395131   50.3237260605 43.4247124698
> > -602501 -1034001
> > 16.3561410682   50.3237171261 43.4247234957
> > -602500 -1034000
> > 16.3561534686   50.3237270561 43.4246990113
> > -602500 -1034001
> > 16.3561550237   50.3237181217 43.4247100363
> >
> >
> > ale gdalinfo dava ulpne jine vysledky:
> > $ gdalinfo pokus_1034001-1034000_602501-602500.tif
> > Driver: GTiff/GeoTIFF
> > Files: pokus_1034001-1034000_602501-602500.tif
> > Size is 2, 2
> > Coordinate System is:
> > PROJCS["unnamed",
> >    GEOGCS["unnamed ellipse",
> >        DATUM["unknown",
> >            SPHEROID["unnamed",6378137,0]],
> >        PRIMEM["Greenwich",0],
> >        UNIT["degree",0.0174532925199433]],
> >    PROJECTION["Mercator_1SP"],
> >    PARAMETER["central_meridian",0],
> >    PARAMETER["scale_factor",1],
> >    PARAMETER["false_easting",0],
> >    PARAMETER["false_northing",0],
> >    UNIT["metre",1,
> >        AUTHORITY["EPSG","9001"]]]
> > Origin = (1820911.314287560060620,6502659.680196234025061)
> > Pixel Size = (1.565026523937731,-1.565026523937731)
> > Metadata:
> >  AREA_OR_POINT=Area
> > Image Structure Metadata:
> >  INTERLEAVE=BAND
> > Corner Coordinates:
> > Upper Left  ( 1820911.314, 6502659.680) ( 16d21'27.09"E, 50d19'28.10"N)
> > Lower Left  ( 1820911.314, 6502656.550) ( 16d21'27.09"E, 50d19'28.04"N)
> > Upper Right ( 1820914.444, 6502659.680) ( 16d21'27.19"E, 50d19'28.10"N)
> > Lower Right ( 1820914.444, 6502656.550) ( 16d21'27.19"E, 50d19'28.04"N)
> > Center      ( 1820912.879, 6502658.115) ( 16d21'27.14"E, 50d19'28.07"N)
> >
> > To jest:
> > Upper Left:  16.3575246439  50.3244732648
> > Lower Left:  16.3575246439  50.3244553137
> > Upper Right: 16.3575527612  50.3244732648
> > Lower Right: 16.3575527612  50.3244553137
> >
> > Pri takhle malem vyrezu by mela byt cisla prakticky stejna, ale ono je to
> > uplne mimo.
> >
> > Michal
> > _______________________________________________
> > FreeGeoCZ mailing list
> > FreeGeoCZ AT fsv.cvut.cz
> > http://mailman.fsv.cvut.cz/mailman/listinfo/freegeocz