EOTR háló készítése

Rövid elkészítési terv:

  1. Paraméterezett rács készítése terjedelemre
  2. Azon téglalapok megtartása, amelyek metszik Magyarországot
  3. Rácsháló téglalapjaihoz tartozó középpontok szelvényszámának megállapítása

QGIS-ben készítsünk rácshálót a következő kóddal (EOV vetületben). A készítendő hálók terjedeleme a következő lesz (xmin, xmax, ymin, ymax): (384000, 960000, 32000, 384000)

Először írjunk egy a kódot a 1:100 000-es méretarányú szelvények elkészítésére. Ekkor a szelvények vízszintes felbontása 48 000, a függőleges pedig 32 000 méter:

from qgis.core import *
import processing

path = "d:/indikatrix-org/eotr-halo-keszitese/" 

filename = "eotr_100000"

extent = "384000, 960000, 32000, 384000"

#processing.algorithmHelp("qgis:creategrid")
parameters = {'TYPE': 2,
              'EXTENT': extent,
              'HSPACING': 48000,
              'VSPACING': 32000,
              'HOVERLAY': 0,
              'VOVERLAY': 0,
              'CRS': 'EPSG:23700',
              'OUTPUT': path + filename + ".shp"
              }
processing.run("qgis:creategrid", parameters)

vectorlayer = QgsVectorLayer(path + filename + ".shp", filename, "ogr")
QgsProject.instance().addMapLayer(vectorlayer)

Fejlesszük tovább a kódot a magyarországi metszet elkészítéséhez:

filename2 = "adm_magyarorszag_eov"
hun_layer = QgsVectorLayer(path + filename2 + ".shp", filename2, "ogr")
#processing.algorithmHelp("native:selectbylocation")
parameters = {'INPUT': gridlayer,
              'PREDICATE': 0,
              'INTERSECT': hun_layer,
              'METHOD': 0
              }
processing.run("native:selectbylocation", parameters)

1. ábra: Szelektált 1:100 000-es EOTR szelvények Magyarországon.
1. ábra: Szelektált 1:100 000-es EOTR szelvények Magyarországon.

Újra fejlesszük tovább kódot, mentsük ki a szelektált elemeket egy új rétegbe:

#processing.algorithmHelp("native:saveselectedfeatures")
parameters = {'INPUT': gridlayer,
              'OUTPUT': path + filename + "_hun" + ".shp"
              }
processing.run("native:saveselectedfeatures", parameters)

gridlayer_hun = QgsVectorLayer(path + filename + "_hun" + ".shp", filename, "ogr")
QgsProject.instance().addMapLayer(gridlayer_hun)

2. ábra: Új rétegbe mentett 1:100 000-es EOTR szelvények Magyarországon.
2. ábra: Új rétegbe mentett 1:100 000-es EOTR szelvények Magyarországon.

Készítsünk egy kódot, ami kiszámítja egy adott ponthoz adott méretarányhoz tartozó EOTR szelvényszámot:

eotrlevels = {
  100000 : {
    "depth" : 0,
    "y" : 48000,
    "x" : 32000
  },
  50000 : {
    "depth" : 1,
    "y" : 24000,
    "x" : 16000
  },
  25000 : {
    "depth" : 2,
    "y" : 12000,
    "x" : 8000
  },
  10000 : {
    "depth" : 3,
    "y" : 6000,
    "x" : 4000
  },
  4000 : {
    "depth" : 4,
    "y" : 3000,
    "x" : 2000
  },
  2000 : {
    "depth" : 5,
    "y" : 1500,
    "x" : 1000
  },
  1000 : {
    "depth" : 6,
    "y" : 750,
    "x" : 500
  }
} 

def quadtree(yf, xf, depth, section_id):
    for i in range(1, depth + 1):
        if i % 3 == 1:
            section_id += "-"        
        c = 2**-i
        if yf >= c:
            # section 2 or 4
            if xf >= c:
                section_id += "2"
            else:
                section_id += "4"
        else:
            # section 1 or 3
            if xf >= c:
                section_id += "1"
            else:
                section_id += "3"
        if yf >= c:
            yf -= c
        if xf >= c:
            xf -= c
    #print(section_id)
    return str(section_id)

def eotr(y, x, m): 
    yn, yr  = divmod(y - 384000, 48000)
    xn, xr  = divmod(x - 32000, 32000)
    depth = eotrlevels[m]['depth']
    section_id = str(int(xn)) + str(int(yn))
    if depth == 0:
        #print(section_id)
        return section_id
    else:
        yf = yr / 48000
        xf = xr / 32000
        section_id = quadtree(yf, xf, depth, section_id)
        return section_id

A kódot a következőképp tudjuk meghívni, tesztelni (pl. Tihanyra):

eotr(562001, 174496, 1000)

Eredmény: 43-412-234

Illesszük ezt a kódot a korábbi kódunk elé, hogy betöltődjenek a függvények, végezzünk apró módosításokat a kódon, majd futtassuk:

from qgis.core import *
from PyQt5.QtCore import *
import processing

eotrlevels = {
  100000 : {
    "depth" : 0,
    "y" : 48000,
    "x" : 32000
  },
  50000 : {
    "depth" : 1,
    "y" : 24000,
    "x" : 16000
  },
  25000 : {
    "depth" : 2,
    "y" : 12000,
    "x" : 8000
  },
  10000 : {
    "depth" : 3,
    "y" : 6000,
    "x" : 4000
  },
  4000 : {
    "depth" : 4,
    "y" : 3000,
    "x" : 2000
  },
  2000 : {
    "depth" : 5,
    "y" : 1500,
    "x" : 1000
  },
  1000 : {
    "depth" : 6,
    "y" : 750,
    "x" : 500
  }
} 

def quadtree(yf, xf, depth, section_id):
    for i in range(1, depth + 1):
        if i % 3 == 1:
            section_id += "-"        
        c = 2**-i
        if yf >= c:
            # section 2 or 4
            if xf >= c:
                section_id += "2"
            else:
                section_id += "4"
        else:
            # section 1 or 3
            if xf >= c:
                section_id += "1"
            else:
                section_id += "3"
        if yf >= c:
            yf -= c
        if xf >= c:
            xf -= c
    #print(section_id)
    return str(section_id)

def eotr(y, x, m): 
    yn, yr  = divmod(y - 384000, 48000)
    xn, xr  = divmod(x - 32000, 32000)
    depth = eotrlevels[m]['depth']
    section_id = str(int(xn)) + str(int(yn))
    if depth == 0:
        #print(section_id)
        return section_id
    else:
        yf = yr / 48000
        xf = xr / 32000
        section_id = quadtree(yf, xf, depth, section_id)
        return section_id
    
#eotr(936000.0, 272000.0, 50000)

path = "d:/indikatrix-org/eotr-halo-keszitese/" 

filename = "eotr_100000"

#processing.algorithmHelp("qgis:creategrid")
extent = "384000, 960000, 32000, 384000"

#processing.algorithmHelp("qgis:creategrid")
parameters = {'TYPE': 2,
              'EXTENT': extent,
              'HSPACING': 48000,
              'VSPACING': 32000,
              'HOVERLAY': 0,
              'VOVERLAY': 0,
              'CRS': 'EPSG:23700',
              'OUTPUT': path + filename + ".shp"
              }
processing.run("qgis:creategrid", parameters)

gridlayer = QgsVectorLayer(path + filename + ".shp", filename, "ogr")
#QgsProject.instance().addMapLayer(gridlayer)

filename2 = "adm_magyarorszag_eov"
hun_layer = QgsVectorLayer(path + filename2 + ".shp", filename2, "ogr")
#processing.algorithmHelp("native:selectbylocation")
parameters = {'INPUT': gridlayer,
              'PREDICATE': 0,
              'INTERSECT': hun_layer,
              'METHOD': 0
              }
processing.run("native:selectbylocation", parameters)
#QgsProject.instance().addMapLayer(gridlayer)

#processing.algorithmHelp("native:saveselectedfeatures")
parameters = {'INPUT': gridlayer,
              'OUTPUT': path + filename + "_hun" + ".shp"
              }
processing.run("native:saveselectedfeatures", parameters)

gridlayer_hun = QgsVectorLayer(path + filename + "_hun" + ".shp", filename, "ogr")

newField = QgsField('label', QVariant.String)
gridlayer_hun.dataProvider().addAttributes([newField])
gridlayer_hun.updateFields()
idx = gridlayer_hun.fields().lookupField('label')

sections = [feature for feature in gridlayer_hun.getFeatures()]
for section in sections:
    y_eov = section.geometry().centroid().asPoint().x()
    x_eov = section.geometry().centroid().asPoint().y()
    eotr_label = eotr(y_eov, x_eov, 100000)
    attrs = { idx : eotr_label }
    fid = section.id()
    gridlayer_hun.dataProvider().changeAttributeValues({ fid : attrs })

QgsProject.instance().addMapLayer(gridlayer_hun)

Végül töltsük be és címkézzük az eredményt a 'label' mező alapján:

3. ábra: 87 darab 1:100 000-es, 48 x 32 km-es EOTR szelvény.
3. ábra: 87 darab 1:100 000-es, 48 x 32 km-es EOTR szelvény.

Módosítsuk úgy a kódot, hogy valamennyi méretarányon elkészüljenek a szelvények, hozzájuk tartozó azonosítókkal (a kód futási ideje a 20 - 30 percet is elérheti):

from qgis.core import *
from PyQt5.QtCore import *
import processing

eotrlevels = {
  100000 : {
    "depth" : 0,
    "y" : 48000,
    "x" : 32000
  },
  50000 : {
    "depth" : 1,
    "y" : 24000,
    "x" : 16000
  },
  25000 : {
    "depth" : 2,
    "y" : 12000,
    "x" : 8000
  },
  10000 : {
    "depth" : 3,
    "y" : 6000,
    "x" : 4000
  },
  4000 : {
    "depth" : 4,
    "y" : 3000,
    "x" : 2000
  },
  2000 : {
    "depth" : 5,
    "y" : 1500,
    "x" : 1000
  },
  1000 : {
    "depth" : 6,
    "y" : 750,
    "x" : 500
  }
} 

def quadtree(yf, xf, depth, section_id):
    for i in range(1, depth + 1):
        if i % 3 == 1:
            section_id += "-"        
        c = 2**-i
        if yf >= c:
            # section 2 or 4
            if xf >= c:
                section_id += "2"
            else:
                section_id += "4"
        else:
            # section 1 or 3
            if xf >= c:
                section_id += "1"
            else:
                section_id += "3"
        if yf >= c:
            yf -= c
        if xf >= c:
            xf -= c
    #print(section_id)
    return str(section_id)

def eotr(y, x, m): 
    yn, yr  = divmod(y - 384000, 48000)
    xn, xr  = divmod(x - 32000, 32000)
    depth = eotrlevels[m]['depth']
    section_id = str(int(xn)) + str(int(yn))
    if depth == 0:
        #print(section_id)
        return section_id
    else:
        yf = yr / 48000
        xf = xr / 32000
        section_id = quadtree(yf, xf, depth, section_id)
        return section_id
    
#eotr(936000.0, 272000.0, 50000)

path = "d:/indikatrix-org/eotr-halo-keszitese/" 

for level in eotrlevels:
    filename = "eotr_" + str(level)
    extent = "384000, 960000, 32000, 384000"

    #processing.algorithmHelp("qgis:creategrid")
    parameters = {'TYPE': 2,
                  'EXTENT': extent,
                  'HSPACING': eotrlevels[level]['y'],
                  'VSPACING': eotrlevels[level]['x'],
                  'HOVERLAY': 0,
                  'VOVERLAY': 0,
                  'CRS': 'EPSG:23700',
                  'OUTPUT': path + filename + ".shp"
                  }
    processing.run("qgis:creategrid", parameters)

    gridlayer = QgsVectorLayer(path + filename + ".shp", filename, "ogr")
    #QgsProject.instance().addMapLayer(gridlayer)

    filename2 = "adm_magyarorszag_eov"
    hun_layer = QgsVectorLayer(path + filename2 + ".shp", filename2, "ogr")
    #processing.algorithmHelp("native:selectbylocation")
    parameters = {'INPUT': gridlayer,
                  'PREDICATE': 0,
                  'INTERSECT': hun_layer,
                  'METHOD': 0
                  }
    processing.run("native:selectbylocation", parameters)
    #QgsProject.instance().addMapLayer(gridlayer)

    #processing.algorithmHelp("native:saveselectedfeatures")
    parameters = {'INPUT': gridlayer,
                  'OUTPUT': path + filename + "_hun" + ".shp"
                  }
    processing.run("native:saveselectedfeatures", parameters)

    gridlayer_hun = QgsVectorLayer(path + filename + "_hun" + ".shp", filename, "ogr")

    newField = QgsField('label', QVariant.String)
    gridlayer_hun.dataProvider().addAttributes([newField])
    gridlayer_hun.updateFields()
    idx = gridlayer_hun.fields().lookupField('label')

    sections = [feature for feature in gridlayer_hun.getFeatures()]
    for section in sections:
        y_eov = section.geometry().centroid().asPoint().x()
        x_eov = section.geometry().centroid().asPoint().y()
        eotr_label = eotr(y_eov, x_eov, level)
        attrs = { idx : eotr_label }
        fid = section.id()
        gridlayer_hun.dataProvider().changeAttributeValues({ fid : attrs })

    QgsProject.instance().addMapLayer(gridlayer_hun)

Ezt követően megállapítható, hogy 249 979 darab 1:1000-es szelvény található Magyarországon. Ha nem tudnánk Magyarország területét, akkor a szelvények méretéből és számából - mint alternatív területszámítási módszerrel - egész jó becslést adhatnánk arra.

4. ábra: Egyre kisebb méretű EOTR szelvények Magyarországon.
4. ábra: Egyre kisebb méretű EOTR szelvények Magyarországon.

5. ábra: EOTR számozás alakulása a 03-as azonosítójú 100 000-es szelvényen, 1:4000 méretarányig, EOV vetület.
5. ábra: EOTR számozás alakulása a 03-as azonosítójú 100 000-es szelvényen, 1:4000 méretarányig.

6. ábra: 06-os azonosítójú szelvény, főképp Szerbiai területtel - Zombor város középen, OpenStreetMap alaptérképpel, Web Mercator vetület.
6. ábra: 06-os azonosítójú szelvény, főképp Szerbiai területtel - Zombor város középen, OpenStreetMap alaptérképpel, Web Mercator vetület.

© Dr. Wirth Ervin