Konfigurace PostGISu pro souřadný systém S-JTSK

V poslední době si tak trochu hraju (no, vlastně trochu víc, je na tom tak nějak postavené téma mé diplomové práce, ale uznávám, že jsem přinejlepším „poučený laik“) s rozšířením RDBMS PostgreSQL PostGIS.

Kromě toho, že se stále přesvědčuju, že odhodit MySQL a plně přejít na Postgres byla dobrá volba (což je vcelku zajímavé, ono už po měsící, kdy používáte namísto MySQL nějakou normální databázi, tak nějak není cesty zpět), musím také semtam řešit nějaké ty problémy.

Zdrojem problémů z nejoblíbenějších jsou převody mezi souřadnými systémy. Pokud už získáte nějaká data pro ČR, tak souřadnice bývají v S-JTSK. Na druhou stranu například pro komunikaci s různými webovými mapami je zapotřebí používat světový standard WGS84. V podstatě do dnešní doby jsme používali poměrně komplikovanou cestu: Nejdříve převedli S-JTSK do S42, potom dále jakýmsi podivným PHP skriptem do WGS84. Naštěstí to nějak fungovalo.

Ale jak jsem říkal, hraju si s PostGISem. A ten jako správny GIS umí (prostřednictvím projektu proj) provádět převody mezi souřadnými systémy, které má popsané.

I přesto, že jsem v tabulce spatial_ref_sys (kde jsou definice jednotlivých souřadných systémů uloženy) nalezl několik systémů, obsahujících zkratku S-JTSK, použití ani jednoho nedávalo uspokojivé výsledky (některé nefungovaly vůbec, jiné se nestrefovaly o pár tisíc kilometrů).

No, není nad to trošku pohledat, a také najít. Bohužel ani použití odkazované definice nepřineslo naprosto dokonalé výsledky (tedy pokud se nesmíříte s chybou v řádu desítek metrů). Chtělo to další hledání, tentokráte poněkud časově náročnější, protože vysněný zdroj se ukrýval pod koncovkou .ppt. Ještě že už jsou dnes vyhledávače tak daleko, že indexují i binární paskvil. Tato práce uvádí úpravy, které zvyšují přesnost převodu. Neměl jsem šanci na nějaké zcela exaktní ověření, nicméně primitivní postup „vzít lokaci adresy v S-JTSK, převést ji a podívat se, jestli na google maps bude opravdu šipka na požadovaném objektu“ poukázal na dostatečnou přesnost převodu.

Definice systému na základě dvou uvedených zdrojů tedy vypadá třeba takto (do pole proj4text jsem za čárky musel vložit mezery, jinak zápis rozhodil zobrazení webu):

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) VALUES (102067, 'ESRI', 102067, 'PROJCS["S-JTSK_Krovak_East_North",GEOGCS["GCS_S_JTSK", DATUM["D_S_JTSK",SPHEROID["Bessel_1841",6377397.155,299.1528128]], PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]], PROJECTION["Krovak"],PARAMETER["False_Easting",0], PARAMETER["False_Northing",0], PARAMETER["Pseudo_Standard_Parallel_1",78.5], PARAMETER["Scale_Factor",0.9999], PARAMETER["Azimuth",30.2881397222222], PARAMETER["Longitude_Of_Center",24.83333333333333], ARAMETER["Latitude_Of_Center",49.5], PARAMETER["X_Scale",-1],PARAMETER["Y_Scale",1], PARAMETER["XY_Plane_Rotation",90],UNIT["Meter",1]]', '+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.2881397222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +towgs84=570.83789, 85.682641,462.84673, 4.9984501,1.5867074,5.2611106,3.5610256');

Sám ji používám ve spojení s databází UIR-ADR, poskytovanou českou státní správou, která obsahuje souřadnice pro krajská města a bývalá okresní města středočeského kraje (a slibují více). Jen je nutné dát pozor, že při definici třeba přes WKT je nutné nejdříve zadat souřadnici y, poté X, obě navíc jako záporné hodnoty. Můj postup pro vytvoření sloupce location v tabulce adresa, který obsahuje souřadnice ve WGS84, je následovný (omlouvám se za případná prasectva a předem díky za opravy):

SELECT AddGeometryColumn('adresa', 'location', 4326, 'POINT', 2);
UPDATE adresa SET location = Transform(GeomFromText('POINT(-' || y || ' -' || x || ')', 102067), 4326);