2

Context

I'd like to download North Rhine-Westphalia's cadastre parcels as GeoJSON. Their website provides a WFS URL:

https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS

I can connect to this endpoint via QGIS. I select the Flurstück (German for parcel) layer and add it:

enter image description here

Result:

enter image description here

Parcels visualized. I'd like to have a local copy of this parcel dataset.


My attempts

Attempt 1: QGIS

I add OpenStreetMap as an XYZ layer to QGIS canvas.

Right click on Flurstück, Export, Save Features As.... I give a filename, choose CRS to be Default CRS: EPSG:4326 - WGS 84, click Extent checkbox, then Calculate from Layer, OpenStreetMap to make sure everything is included. I set coordinate precision to 8.

enter image description here

Explore the output using GeoPandas:

import geopandas as gpd
df = gpd.read_file("NRW_parcels.geojson")

print(len(df)) # outputs: 100000

print(df.geometry)

The geometries exported are Polygons, as expected:

0        POLYGON ((6.39946 50.34495, 6.39977 50.34465, ...
1        POLYGON ((6.40553 50.32334, 6.40559 50.32333, ...
2        POLYGON ((6.40965 50.33726, 6.41024 50.33669, ...
3        POLYGON ((6.40781 50.33902, 6.40820 50.33865, ...
4        POLYGON ((6.40380 50.34110, 6.40615 50.34167, ...
                               ...                        
99995    POLYGON ((6.85390 50.53447, 6.85393 50.53442, ...
99996    POLYGON ((6.88087 50.53661, 6.88093 50.53629, ...
99997    POLYGON ((6.85023 50.53525, 6.85109 50.53488, ...
99998    POLYGON ((6.85343 50.53409, 6.85356 50.53364, ...
99999    POLYGON ((6.88043 50.53765, 6.88049 50.53743, ...
Name: geometry, Length: 100000, dtype: geometry

However, there are only 100000 of them. There are more parcels than that in NRW, so we failed to export all of them.


Attempt 2: ogr2ogr with URL

Load the WFS URL into a variable:

URL="https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS"

Explore options with ogrinfo:

ogrinfo WFS:$URL

Response:

INFO: Open of `WFS:https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS'
      using driver `WFS' successful.
Metadata:
  ABSTRACT=Dieser Dienst stellt für das INSPIRE-Thema Flurstücke/Grundstücke (Katasterparzellen) aus ALKIS umgesetzte Daten bereit.
  PROVIDER_NAME=Geobasis NRW
  TITLE=INSPIRE-WFS NW Flurstücke/Grundstücke ALKIS
ERROR 1: No schema information loaded
1: cp:CadastralParcel (title: Flurstück) (Point)
ERROR 1: No schema information loaded
2: cp:CadastralZoning (title: Katasterbezirk) (Point)

Explore the cp:CadastralParcel layer specifically:

ogrinfo WFS:$URL cp:CadastralParcel

Long output. Beginning of it:

INFO: Open of `WFS:https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS'
      using driver `WFS' successful.
Metadata:
  ABSTRACT=Dieser Dienst stellt für das INSPIRE-Thema Flurstücke/Grundstücke (Katasterparzellen) aus ALKIS umgesetzte Daten bereit.
  PROVIDER_NAME=Geobasis NRW
  TITLE=INSPIRE-WFS NW Flurstücke/Grundstücke ALKIS
ERROR 1: No schema information loaded

Layer name: cp:CadastralParcel
Metadata:
  ABSTRACT=Gebiete, die anhand des Liegenschaftskatasters oder gleichwertiger Verzeichnisse bestimmt werden.
  KEYWORD_1=Flurstück
  KEYWORD_2=Cadastral Parcel
  TITLE=Flurstück
Geometry: Point
Feature Count: 9355266
ERROR 1: No schema information loaded
Extent: (300140.721000, 5577691.422000) - (353253.838000, 5605660.731000)
Layer SRS WKT:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
        BBOX[38.76,6,84.33,12.01]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
gml_id: String (0.0)
identifier: String (0.0)
areaValue: Real (0.0)
beginLifespanVersion: String (0.0)
endLifespanVersion: String (0.0)
localId: String (0.0)
namespace: String (0.0)
label: Integer (0.0)
nationalCadastralReference: String (0.0)
validFrom: String (0.0)
validTo: String (0.0)
zoning: String (0.0)
OGRFeature(cp:CadastralParcel):1
  gml_id (String) = CadastralParcel_05439400700069______
  identifier (String) = https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400700069______
  areaValue (Real) = 35190
  beginLifespanVersion (String) = 2017-01-16T15:22:59Z
  endLifespanVersion (String) = 
  localId (String) = CadastralParcel_05439400700069______
  namespace (String) = https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis
  label (Integer) = 69
  nationalCadastralReference (String) = 05439400700069______
  validFrom (String) = 1974-01-01
  validTo (String) = 
  zoning (String) = 
  POINT (315035.325 5579575.262)

OGRFeature(cp:CadastralParcel):2
  gml_id (String) = CadastralParcel_05439400800001______
  identifier (String) = https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400800001______
  areaValue (Real) = 21530
  beginLifespanVersion (String) = 2021-04-08T11:29:11Z
  endLifespanVersion (String) = 
  localId (String) = CadastralParcel_05439400800001______
  namespace (String) = https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis
  label (Integer) = 1
  nationalCadastralReference (String) = 05439400800001______
  validFrom (String) = 1974-01-01
  validTo (String) = 
  zoning (String) = 
  POINT (315428.343 5578357.297)

We learnt: there are 9355266 parcels.

Attempt to download them:

rm NRW_parcels.geojson
ogr2ogr -f GeoJSON NRW_parcels.geojson WFS:$URL 

Output:

ERROR 1: No schema information loaded
ERROR 1: No schema information loaded
Warning 1: Value '54/2' of field cp:CadastralParcel.label parsed incompletely to integer 54.
Warning 1: Value '33/9' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/7' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/12' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/6' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/11' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/8' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '33/10' of field cp:CadastralParcel.label parsed incompletely to integer 33.
Warning 1: Value '3/2' of field cp:CadastralParcel.label parsed incompletely to integer 3.
Warning 1: Value '54/1' of field cp:CadastralParcel.label parsed incompletely to integer 54.
Warning 1: Value '3/1' of field cp:CadastralParcel.label parsed incompletely to integer 3.
Warning 1: Value '48/1' of field cp:CadastralParcel.label parsed incompletely to integer 48.
Warning 1: Value '48/2' of field cp:CadastralParcel.label parsed incompletely to integer 48.
Warning 1: Value '37/1' of field cp:CadastralParcel.label parsed incompletely to integer 37.

... # warnings go on for long

More than 1000 errors or warnings have been reported. No more will be reported from now.

This produces the ~60 MB NRW_parcels.geojson geojson file. Check how it looks like:

head -c2000 NRW_parcels.geojson

Output:

{
"type": "FeatureCollection",
"name": "cp:CadastralParcel",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::25832" } },
"features": [
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400700069______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400700069______", "areaValue": 35190.0, "beginLifespanVersion": "2017-01-16T15:22:59Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400700069______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 69, "nationalCadastralReference": "05439400700069______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordinates": [ 315035.325000000011642, 5579575.262000000104308 ] } },
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400800001______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400800001______", "areaValue": 21530.0, "beginLifespanVersion": "2021-04-08T11:29:11Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400800001______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 1, "nationalCadastralReference": "05439400800001______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordinates": [ 315428.343, 5578357.29700000025332 ] } },
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400700065______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400700065______", "areaValue": 2702.0, "beginLifespanVersion": "2017-01-16T15:22:59Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400700065______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 65, "nationalCadastralReference": "05439400700065______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordin

Unfortunately, these are just points, not Polygons. Try specify the layer we want:

rm NRW_parcels.geojson
ogr2ogr -f GeoJSON NRW_parcels.geojson WFS:$URL cp:CadastralParcel

Same result as before.

Check how many points we have via Python:

df = gpd.read_file("NRW_parcels.geojson")
print(len(df)) # outputs: 100000

Problems: instead of Polygons we have Points, instead of 9355266 of them we have only 100000. Try fixing this latter problem using this answer to the question How to make a local copy of a large WFS dataset? :

rm NRW_parcels.geojson
ogr2ogr -f GeoJSON NRW_parcels.geojson -dsco OGR_WFS_PAGING_ALLOWED=ON WFS:$URL cp:CadastralParcel 

Output:

Warning 6: driver GeoJSON does not support creation option OGR_WFS_PAGING_ALLOWED

then the same as before.


Attempt 3: ogr2ogr with xml

This answer to the question How to make a local copy of a large WFS dataset? describes how to use an xml file with ogr2ogr to access data via a WFS server. Create NRW_parcels.xml with these contents:

<OGRWFSDataSource>
  <URL>https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS</URL>
  <PagingAllowed>ON</PagingAllowed>
</OGRWFSDataSource>  

Attempt ogr2ogr:

rm NRW_parcels.geojson
ogr2ogr -f GeoJSON NRW_parcels.geojson NRW_parcels.xml

Same behaviour as before.


Attempt 4: wget

This answer to the question How to make a local copy of a large WFS dataset? suggests wget:

rm NRW_parcels.geojson
wget -O NRW_parcels.geojson $URL

Output:

--2022-11-30 11:04:50--  https://www.wfs.nrw.de/geobasis/wfs_nw_inspire-flurstuecke_alkis?REQUEST=GetCapabilities&SERVICE=WFS
Resolving www.wfs.nrw.de (www.wfs.nrw.de)... 93.184.133.187
Connecting to www.wfs.nrw.de (www.wfs.nrw.de)|93.184.133.187|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: unspecified [text/xml]
Saving to: ‘NRW_parcels.geojson’

NRW_parcels.geojson                     [ <=>                                                               ]  22,60K  --.-KB/s    in 0,05s   

2022-11-30 11:04:50 (454 KB/s) - ‘NRW_parcels.geojson’ saved [23143]

A 23 KB NRW_parcels.geojson was created, top part via head -c2000 NRW_parcels.geojson:

<?xml version="1.0" encoding="utf-8"?>
<wfs:WFS_Capabilities xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:wfs="http://www.opengis.net/wfs/2.0" xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:fes="http://www.opengis.net/fes/2.0" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:inspire_common="http://inspire.ec.europa.eu/schemas/common/1.0" xmlns:inspire_dls="http://inspire.ec.europa.eu/schemas/inspire_dls/1.0" xsi:schemaLocation="http://www.opengis.net/wfs/2.0 http://www.wfs.nrw.de/aaa-suite/schema/ogc/wfs/2.0/wfs.xsd http://inspire.ec.europa.eu/schemas/inspire_dls/1.0 http://www.wfs.nrw.de/aaa-suite/schema/inspire/inspire_dls/1.0/inspire_dls.xsd" version="2.0.0">
<ows:ServiceIdentification>
<ows:Title>INSPIRE-WFS NW Flurstücke/Grundstücke ALKIS</ows:Title>
<ows:Abstract>Dieser Dienst stellt für das INSPIRE-Thema Flurstücke/Grundstücke (Katasterparzellen) aus ALKIS umgesetzte Daten bereit.</ows:Abstract>
<ows:Keywords>
<ows:Keyword>Geobasisdaten</ows:Keyword>
<ows:Keyword>AdV</ows:Keyword>
<ows:Keyword>WFS</ows:Keyword>
<ows:Keyword>AdV-OWS-Basisprofil</ows:Keyword>
<ows:Keyword>AdV-WFS-Profil 2.0.0</ows:Keyword>
<ows:Keyword>AdV-INSPIRE-Produktspezifikation (Version 1.0.0)</ows:Keyword>
<ows:Keyword>NW</ows:Keyword>
<ows:Keyword>Nordrhein-Westfalen</ows:Keyword>
<ows:Keyword>inspireidentifiziert</ows:Keyword>
<ows:Keyword>infoFeatureAccessService</ows:Keyword>
<ows:Keyword>INSPIRE-WFS_NW_Flurstuecke_Grundstuecke_ALKIS</ows:Keyword>
<ows:Keyword>ALKIS</ows:Keyword>
<ows:Keyword>Amtliches Liegenschaftskatasterinformationssystem</ows:Keyword>
<ows:Keyword>Flurstücke/Grundstücke</ows:Keyword>
<ows:Keyword>Cadastral Parcels</ows:Keyword>
</ows:Keywords>
<ows:ServiceType>WFS</ows:ServiceType>
<ows:ServiceTypeVersion>2.0.0</ows:ServiceTypeVersion>
<ows:Fees>Nutzungsbedingungen: Es gelten die durch den IT-Planungsrat im Datenportal f&#252;r Deutschland (GovData) ve%                        

Fail.


Attempt 5: ogr2ogr with paging

This answer to the question Save WFS (~100000 features) as geopackage with QGIS uses startindex and count in the URL.

Try:

rm NRW_parcels.geojson 
ogr2ogr -f GeoJSON NRW_parcels.geojson "WFS:${URL}&startindex=1000&count=2000" cp:CadastralParcel

Output:

ERROR 1: No schema information loaded
ERROR 1: No schema information loaded

A 1.2 MB NRW_parcels.geojson is produced. Check the beginning of the file via head -c2000 NRW_parcels.geojson:

{
"type": "FeatureCollection",
"name": "cp:CadastralParcel",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::25832" } },
"features": [
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400700069______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400700069______", "areaValue": 35190.0, "beginLifespanVersion": "2017-01-16T15:22:59Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400700069______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 69, "nationalCadastralReference": "05439400700069______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordinates": [ 315035.325000000011642, 5579575.262000000104308 ] } },
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400800001______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400800001______", "areaValue": 21530.0, "beginLifespanVersion": "2021-04-08T11:29:11Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400800001______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 1, "nationalCadastralReference": "05439400800001______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordinates": [ 315428.343, 5578357.29700000025332 ] } },
{ "type": "Feature", "properties": { "gml_id": "CadastralParcel_05439400700065______", "identifier": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis/CadastralParcel_05439400700065______", "areaValue": 2702.0, "beginLifespanVersion": "2017-01-16T15:22:59Z", "endLifespanVersion": "", "localId": "CadastralParcel_05439400700065______", "namespace": "https://registry.gdi-de.org/id/de.nw.inspire.cp.alkis", "label": 65, "nationalCadastralReference": "05439400700065______", "validFrom": "1974-01-01", "validTo": "", "zoning": "" }, "geometry": { "type": "Point", "coordin%

Still only Points unfortunately. Check their number via Python:

import geopandas as gpd
df = gpd.read_file("NRW_parcels.geojson")
print(len(df)) # outputs: 2000

Probably I could loop through all of them then, by increasing startindex, and get many Points. Those are still not Polygons though, and I'd like to have Polygons.


Question

How can I download all 9355266 parcel Polygons of North Rhine-Westphalia?

8
  • 1
    When I load your WFS layer in QGIS, even though the "feature count" shows 9M entries, the Attribute Table shows exactly 100000 entries. It looks like there is a limit on the WFS server, which returns an error " Parameter 100000 exceeds configurated limit" if you call it this way in a browser: wfs.nrw.de/geobasis/…. One approach might be do download/merge the data from multiple calls to the WFS server on different bounding boxes
    – Kasper
    Commented Nov 30, 2022 at 11:42
  • 1
    That server only supports GML output - <ows:Parameter name="outputFormat"> <ows:AllowedValues> <ows:Value>text/xml; subtype=gml/3.2.1</ows:Value> <ows:Value>application/gml+xml; version=3.2</ows:Value> </ows:AllowedValues>
    – Ian Turton
    Commented Nov 30, 2022 at 11:49
  • 1
    This task verges on web scraping, which is usually a TOS violation on the web service. I suggest you contact the data host and try to obtain the data directly.
    – Vince
    Commented Nov 30, 2022 at 13:21
  • 2
    Of course GDAL supports paging with page_size gdal.org/drivers/vector/wfs.html. However, you are reading data that follows the INSPIRE data specification for cadastral parcels. That is using an application schema and contains complex features as appears from the data model inspire.ec.europa.eu/data-model/approved/r4618-ir/html/…. You cannot save such data into GeoJSON. I would try to save all the data from the WFS into local GML file and continue then with the GDAL GMLAS driver gdal.org/drivers/vector/gmlas.html.
    – user30184
    Commented Nov 30, 2022 at 18:42
  • 1
    Downloading all data from a WFS is absolutely what a WFS is for, don't ever be afraid about violating TOS by using a service meant for public usage. Commented May 25, 2023 at 12:54

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.