Formule berekening tegelnummer download BGT

[rboeters]

Ik zou graag een formule krijgen om op basis van x, y (RD New) een tegelnummer te kunnen berekenen om zo de link van de BGT download via https://www.pdok.nl/nl/producten/pdok-downloads/download-basisregistratie-grootschalige-topografie te kunnen samenstellen.

Dit om downloads te kunnen automatiseren.

Kunnen jullie deze vrijgeven?

[wouter.visscher]

Hoi @rboeters,

Het betreft de z-order curve (Z-order curve - Wikipedia), het grid is met die logica uitgezet. De BBOX betreft: -100000,200000, 412000,712000 (RD). Het is dus een grid van 512000x512000 meter verdeeld in tegels van 2000x2000 meter, op het laagste niveau.

https://downloads.pdok.nl/service/extract.zip?extractname=bgt&extractset=citygml&excludedtypes=plaatsbepalingspunt&history=true&tiles={"layers"%3A[{"aggregateLevel"%3A1%2C"codes"%3A[4094]}]}

De url heeft de query parameter ‘tiles’ welke een encoded stukje json consumeert.

Bijvoorbeeld (decoded): {“layers”:[{“aggregateLevel”:0,“codes”:[16378]}]}

Dit request vraagt dus op het laagste niveau (0) de tegel met de code 16378 op. De parameter codes is een array, je kan dus meerdere tegels op 1 zoomniveau opvragen door bijvoorbeeld mee te geven:

{“layers”:[{“aggregateLevel”:0,“codes”:[16376, 16377, 16378, 16379]}]}

Maar in het bovenste geval vormen de viertegel op aggregateLevel 1 (een niveau hoger) samen een nieuwe tegel namelijk tegel: 4094. Dus het volgende request levert hetzelfde resultaat op:

{“layers”:[{“aggregateLevel”:1,“codes”:[4094]}]}

Want:
16376/4 → floor(4094) = 4094
16377/4 → floor(4094,25) = 4094
16378/4 → floor(4094,5) = 4094
16379/4 → floor(4094,75) = 4094

Met behulp van deze logica kan je ‘navigeren’ tussen verschillende niveaus, door te ‘delen’ of te ‘vermenigvuldigen’ met 4.

Om trouwens antwoord te geven op je vraag: om van een RD-coord naar een tile number te gaan, dat doen wij door middel van ‘bit shifting’. DEPRECATED_featured-shared/ZOrder.java at master · PDOK/DEPRECATED_featured-shared · GitHub

In Postgresql zou je het volgende kunnen doen trouwens (als dat makkelijker is), dit zou je eenvoudig in een functie kunnen schroeven (ff een disclaimer: code is een voorbeeld en geen ‘officiele’ impementatie, enz, enz…)

DO $$
DECLARE
	morton_id	INTEGER := 0;
	x_masked_i	INTEGER;
	y_masked_i	INTEGER;
	i		INTEGER := 0;
	x_coord		NUMERIC(10,1) := 42459;
	y_coord		NUMERIC(10,1) := 383029;
	x		NUMERIC(10,1);
	y		NUMERIC(10,1);

BEGIN
	x := ((x_coord - -100000)/2000);
	IF ((x % 1) != 0) THEN
		x := FLOOR(x);
	END IF;

	y := ((y_coord - 200000)/2000);
	IF ((y % 1) != 0) THEN
		y := FLOOR(y);
	END IF;		

	WHILE (i < 32) LOOP
		x_masked_i = (x::INTEGER & (1 << i));
		y_masked_i = (y::INTEGER & (1 << i));

		morton_id := morton_id | (x_masked_i << i);
		morton_id := morton_id | (y_masked_i << (i + 1));

		i := i + 1;				
	END LOOP;

	RAISE INFO 'morton_id: %', morton_id;
END
$$;

[emacgillavry]

Dankjewel @wouter.visscher voor dit uitgebreide antwoord, zeg! Bit-shifting is inderdaad de secret sauce in de meeste code rondom tiling schema’s! Leuk om dat hier ook weer terug te zien komen.

[rboeters]

Inderdaad, bedankt @wouter.visscher! Hier gaat het zeker mee lukken.

[Peter_Hh]

Heb de code geschreven in c#. Voor de liefhebbers:

public int tileNumberBGT(double X, double Y)

{
  int tileNr = 0;
  double lowerX = -100000;
  double lowerY = 200000;
  double upperX = 412000;
  double upperY = 712000;
  double tileScale = 2000;

  double x = Math.Floor((X - lowerX) / tileScale);
  double y = Math.Floor((Y - lowerY) / tileScale);
  int iX;
  int iY;
  for (int i = 0; i < 32; i++)
  {
    iX = Convert.ToInt32(x) & (1 << i);
    iY = Convert.ToInt32(y) & (1 << i);
    tileNr = tileNr | (iX << i);
    tileNr = tileNr | (iY << (i + 1));
  }
  return tileNr;
}

[JvdB]

Heeft iemand toevallig ook nog een omgedraaide versie op de planken liggen? Van tilenummer naar minX en minY van de tile :)?

[rboeters]

In C++ heb ik dat… Is niet één op één te gebruiken, maar wellicht heb je wat aan het idee.

CExtents PDOKStufDownloader::getExtentsForTile(int tilenr)
{
  // See https://geoforum.nl/t/formule-berekening-tegelnummer-download-bgt/658 for explanation of 
  // calculation
  double x_lo = 0;
  double y_lo = 0;

  CGGString binary = QString::fromStdString(std::bitset<32>(tilenr).to_string());
  CGGString x_binary = L"";
  CGGString y_binary = L"";

  bool addToX = true;

  for (int i = (binary.length()-1); i > -1; --i)
  {
    if (addToX)
    {
      x_binary = binary.GetAtA(i) + x_binary;
      addToX = false;
    }
    else
    {
      y_binary = binary.GetAtA(i) + y_binary;
      addToX = true;
    }
  }
      
  x_lo = std::stoi(QString(x_binary).toStdWString(), nullptr, 2);
  y_lo = std::stoi(QString(y_binary).toStdWString(), nullptr, 2);

  x_lo = -100000 + (2000*x_lo);
  y_lo = 200000 + (2000*y_lo);

  CExtents extents(x_lo, y_lo, x_lo+2000, y_lo+2000);
  return extents;
}

[JvdB]

Held!

[JvdB]

@rboeters ik kom er toch achter dat ik iets minder in C++ thuis ben dan verwacht haha .

Zou je iets meer kunnen toelichten over de volgende statements doen met bijv tile # 48:

@rboeters

CGGString binary = QString::fromStdString(std::bitset<32>(tilenr).to_string());

Is dit een String byte representatie? tile 48 maakt : “110000”?

@rboeters

binary.GetAtA(i) + x_binary;

Wat doet de getAtA? is dat de waarde van de String representatie op die index?
Dwz met 110000 en getAtA(0) → 1, getAtA(1) → 1, getAtA(2) → 0?

@rboeters

std::stoi(QString(x_binary).toStdWString(), nullptr, 2);

Is dit simpelweg van Stringbyte respresentatie naar double? → “110000” naar 110000?

Alvast bedankt voor het antwoord!

[Edwin]

Hallo allemaal,

Is er iemand die van deze tegels een datasetje heeft, een shapefile :-o of een GeoPackage ?
Bij het downloaden van BGT data vanuit andere software wil je nét niet nog een extra paar tegels selecteren met de polygoon.

Bvd.

[Anton]

InfraCAD Map heeft een XML bestand met alle tegels en de coördinaten.

Maar de tegeltjes gaan in de toekomst verdwijnen. De nieuwe omgeving waar je BGT en BRK gaat downloaden werkt niet meer met tegels.

[wouter.visscher]

Hoi @Edwin,
Hierbij het oude stukje code aangepast dat deze de tiles INSERT in postgis, voor een grid van 256x256.

CREATE TABLE bgt.tiles
(
  id integer,
  geometry geometry(Polygon,28992)
);

DO language plpgsql $$
DECLARE 
  xaxis INTEGER := 0;
  yaxis INTEGER := 0;
  lowx NUMERIC(10,1) := -100000;
  lowy NUMERIC(10,1) :=  200000;
  tilesize NUMERIC(10,1) := 2000;
BEGIN
  WHILE (xaxis< 257) LOOP
    yaxis := 0;
    WHILE (yaxis < 257) LOOP
      DECLARE
        poly TEXT;
        morton_id INTEGER := 0;
        x_coord		NUMERIC(10,1) := 42459;
        y_coord		NUMERIC(10,1) := 383029;			  
      BEGIN

        x_coord := lowx + (xaxis*tilesize) + (tilesize/2) ;
        y_coord := lowy + (yaxis*tilesize) + (tilesize/2);

        DECLARE
          x_masked_i	INTEGER;
          y_masked_i	INTEGER;
          i		INTEGER := 0;
          x		NUMERIC(10,1);
          y		NUMERIC(10,1);			
        BEGIN
          x := ((x_coord - -100000)/2000);
          IF ((x % 1) != 0) THEN
            x := FLOOR            x_coord := lowx + (xaxis*tilesize) + (tilesize/2) ;
        y_coord := lowy + (yaxis*tilesize) + (tilesize/2);rd - 200000)/2000);
          IF ((y % 1) != 0) THEN
            y := FLOOR(y);
          END IF;		

          WHILE (i < 32) LOOP
            x_masked_i = (x::INTEGER & (1 << i));
            y_masked_i = (y::INTEGER & (1 << i));

            morton_id := morton_id | (x_masked_i << i);
            morton_id := morton_id | (y_masked_i << (i + 1));

            i := i + 1;				
          END LOOP;			
        END;
      
        poly := format('POLYGON((%1$s %2$s, %3$s %4$s, %5$s %6$s, %7$s %8$s, %1$s %2$s))', 
            lowx + (xaxis*tilesize), lowy + (yaxis*tilesize), 
            lowx + (xaxis*tilesize) + tilesize, lowy + (yaxis*tilesize), 
            lowx + (xaxis*tilesize) + tilesize, lowy + (yaxis*tilesize) + tilesize, 
            lowx + (xaxis*tilesize), lowy + (yaxis*tilesize) + tilesize);
        insert into bgt.tiles(id, geometry) values (morton_id, ST_GeomFromText(poly, 28992));

      END;
      --RAISE INFO 'xaxis: %, yaxis: %', xaxis, yaxis;
      yaxis := yaxis + 1;				
    END LOOP;
                xaxis := xaxis + 1;	
  END LOOP;
END $$;

En een geopackage van het betreffende grid, voor de liefhebbers.