Extract administrative boundaries from citypopulation.de as GeoJSON

admin.md

Open page, for example, and type the following commands into web browser's console and copy output to file:

script = document.createElement('script');
document.head.appendChild(script);

script.onload = function () {
    polygons = [];
    adminLayer.forEach(function (polygon) {
        linearRing = [];
        polygon.getPath().b.forEach(function (vertex, index) {
            linearRing.push([vertex.lng(), vertex.lat()]);
        });

        firstCoordinate = linearRing[0];
        lastCoordinate = linearRing.slice(-1)[0];
        if ((firstCoordinate[0] !== lastCoordinate[0]) ||
            (firstCoordinate[1] !== lastCoordinate[1])) {
            linearRing.push(firstCoordinate);
        }

        polygons.push(
            new ol.Feature({
                geometry: new ol.geom.Polygon([linearRing]),
                name: polygon.name
            })
        );
    });

    geojson = new ol.format.GeoJSON();
    console.log(geojson.writeFeatures(polygons));
};

script.src = 'https://cdnjs.cloudflare.com/ajax/libs/ol3/4.2.0/ol.js';

I get an error when trying to execute the script:

VM113:6 Uncaught TypeError: adminLayer.forEach is not a function
    at HTMLScriptElement.script.onload (<anonymous>:6:16)

Is the script working now, or have citypopulation.de developers changed something?

We receive the same error, was any of you able to solve the issue?

@jimmylevell it seems people from citypopulation.de want to keep their data in secret. But since we are dealing with javascript, then the decoding algorithm is available. Here is a Python snippet which demonstrates how geometries can be extracted:

import json
from requests import Session
from pyproj import proj
from shapely.geometry import Polygon, MultiPolygon, mapping
from pyproj import Transformer
from pyproj.crs import CRS


country = "luxembourg"

s = Session()
req1 = s.get("https://www.citypopulation.de/en/{}/admin/".format(country))

# mode: ["adm1", "adm2"]
adm = s.get(
    "https://www.citypopulation.de/proc/retrieve_adminareas.php",
    params={
        "reqid": 1,
        "pageid": "{}-admin".format(country),
        "type": "adm1",
        "mode": "all",
        "objid": "",
        "cache": "",
    },
).json()


def decode_coordinates(geom_type, s, precision=None):
    def f(name, m, fn):
        def c(p, a, f, z):
            def g(z, d, t):
                c = 0
                for i in range(len(z)):
                    b = ord(z[i])
                    if b == 40:
                        b = 92
                    if b == 41:
                        b = 63
                    c = c * t + b - d
                return c

            k = None
            if "#" not in a:
                if "$" not in a:
                    if "%" not in a:
                        if "&" not in a:
                            return None
                        i = a.index("&")
                        k = [-1, -1]
                    else:
                        i = a.index("%")
                        k = [-1, 1]
                else:
                    i = a.index("$")
                    k = [1, -1]
            else:
                i = a.index("#")
                k = [1, 1]
            p.append(g(a[:i], f, z) * k[0])
            p.append(g(a[i + 1 :], f, z) * k[1])

        memo = []
        if not isinstance(m, list):
            e = m.split("!")
            m = []
            for i in range(len(e)):
                c(m, e[i], 47, 80)
        word = m[0]
        token = m[1]
        result = by_xy_res(word, token, name, fn)
        memo.append(result)
        for i in range(2, len(m), 2):
            word = word + m[i]
            token = token + m[i + 1]
            result = by_xy_res(word, token, name, fn)
            memo.append(result)
        return memo

    if precision is None:
        precision = 1e5
    coords = []
    geometry = None
    for index in range(len(s)):
        if "-" == s[index]:
            if geom_type == "poly":
                for part in coords:
                    part.append(part[0])
                polygon = Polygon([[p[0], p[1]] for poly in coords for p in poly])
                if geometry is None:
                    geometry = MultiPolygon([polygon])
                else:
                    geometry = MultiPolygon(
                        [geom for geom in geometry.geoms]
                        + [
                            polygon,
                        ]
                    )
                coords = []
        else:
            coords.append(f(precision, s[index], "EPSG:3857"))

    if geom_type == "poly":
        for part in coords:
            part.append(part[0])
        polygon = Polygon([[p[0], p[1]] for poly in coords for p in poly])
        if geometry is None:
            geometry = MultiPolygon([polygon])
        else:
            geometry = MultiPolygon(
                [geom for geom in geometry.geoms]
                + [
                    polygon,
                ]
            )
    return geometry


transformers = {
    "EPSG:3857": Transformer.from_crs(
        CRS.from_epsg(4326), CRS.from_epsg(3857), always_xy=True
    )
}


def by_xy_res(t, n, i, r):
    if i:
        t /= i
        n /= i
    return by_latlon(n, t, r)


def by_latlon(t, n, i):
    transformer = transformers[i]
    if i and i != "EPSG:4326":
        return transformer.transform(n, t)
    else:
        return [n, t]


geoms = []
for feature in adm["objs"]:
    print(feature["id"])
    geoms.append(decode_coordinates("poly", feature["c"]))


geojson = {}
geojson["type"] = "FeatureCollection"
geojson["features"] = [
    {"type": "Feature", "geometry": g} for g in [mapping(p) for p in geoms]
]

with open("/tmp/output.js", "w") as dst:
    json.dump(geojson, dst)

N.B. Although it is technically possible to extract coordinate data of places but the maps and other geospatial data are under the copyright of »City Population« and/or under the copyright of the organizations denoted at the maps (ToU).