User:Twse/ektowld

Lantmäteriet publicerar gamla ekonomiska kartor, de kan man göra tiles av och använda i olika program, se User:Dalkvist/ekonomiska kartan.

I JOSM är ett alternativ att skapa en world file och öppna dem med pluginen PicLayer, enligt instruktion nedan.

Ladda hem karta

1. Ladda ner ett kartblad av Ekonomiska kartan från lantmäteriet

   Varje karta har en beskärd sida och en komplett sida

2. Exportera den beskärda sidan som en PNG-fil
3. På den kompletta sidan står det två nummer längst ner till vänster, notera dem

Skapa world file

1. Ladda ner och installera python och python-gdal
2. Spara skriptet nedan som ektowld.py
3. Kör skriptet med PNG-filen och nummer från kartan som parametrar, ex:

   python ektowld.py karta.png 6730 1565

4. Skriptet skapar en .wld och .prj fil namngivna efter PNG-filen

Öppna i JOSM

1. Aktivera plugin PicLayer
2. Välj PicLayer/New picture layer from file, välj PNG-filen
3. PicLayer kommer automatisk georeferera lagret med hjälp av world-filen

Skript

#!/usr/bin/python
# -*- coding: utf-8 -*-
import os
import sys

import osgeo.gdal as gdal
import osgeo.osr as osr

filnamn = sys.argv[1]

#2x 4 siffror nere till vänster på kartblad
kartKoordY = float(sys.argv[2])*1000
kartKoordX = float(sys.argv[3])*1000

def RT90toWEB(rt90x, rt90y):
    rt90 = osr.SpatialReference()
    #Nyare metod från Lantmäteriet
    rt90.ImportFromProj4("""+proj=tmerc +ellps=GRS80 +a=6378137
                            +rf=298.257222101 +lon_0=15.8062845294
                            +k=1.00000561024 +y_0=-667.711 +x_0=1500064.274
                            +datum=WGS84 +units=m""")
    
    #Web-format använt internt i JOSM/PicLayer
    web = osr.SpatialReference()
    web.ImportFromEPSG(3857)
    
    ct = osr.CoordinateTransformation(rt90, web)
    return ct.TransformPoint(rt90x, rt90y)

def RT90toGCP(rt90x, rt90y, pixelX, pixelY):
    webX, webY, webHeight = RT90toWEB(rt90x, rt90y)
    gcp = gdal.GCP(webX, webY, webHeight, pixelX, pixelY)
    return gcp

gcps=[]

#SydVäst
gcps.append( RT90toGCP(kartKoordX, kartKoordY, 0, 5000)   )
#SydÖst
gcps.append( RT90toGCP(kartKoordX+5000, kartKoordY, 5000, 5000) )
#NordVäst
gcps.append( RT90toGCP(kartKoordX, kartKoordY+5000, 0,0) )
#NordÖst
gcps.append( RT90toGCP(kartKoordX+5000, kartKoordY+5000, 5000, 0) )

xform = gdal.GCPsToGeoTransform(gcps) 

#.wld läses automatiskt av PicLayer
wld = open(os.path.splitext(filnamn)[0] + '.wld', 'w')
wld.write("%0.10f\n" % xform[1])
wld.write("%0.10f\n" % xform[2])
wld.write("%0.10f\n" % xform[4])
wld.write("%0.10f\n" % xform[5])
wld.write("%0.10f\n" % xform[0])
wld.write("%0.10f\n" % xform[3])
wld.close()

#.prj används inte av PicLayer, kan vara bra till andra program
web = osr.SpatialReference()
web.ImportFromEPSG(3857)
webWkt = web.ExportToWkt()
prj = open(os.path.splitext(filnamn)[0] + '.prj', 'w')
prj.write(webWkt)
prj.close()