1

I wish to interrogate a WFS layer: 'https://environment.data.gov.uk/spatialdata/flood-map-for-planning-rivers-and-sea-flood-zone-2/wfs'

I am not the greatest python programmer (not quick) but can do a lot of ESRI related functions although struggle with setting environments, adding new libraries etc, but I am stuck with this WFS problem:

  • I wish to define a bounding box which will change depending on a client site (example here): minx = 419000, miny = 419000, maxx = 421000, maxy = 421000

  • I wish to save the resulting features (polygons) as a shapefile

  • I have all other code in Python 2.7.14 so wish to use same environment

  • I can't find anything online that explains from start to finish

Is anyone out there that can explain / show example of something they have done (different URL is fine) that will point me into the right direction

1
  • Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.
    – Community Bot
    Sep 10, 2021 at 10:02

3 Answers 3

1

You have two steps one download features from WFS and second write a shapefile fortunately both a common operations and are documented in the OGR Cookbook.

So cutting and pasting gives this (You will need to adjust the field names ect):

import sys

try:
    from osgeo import ogr, osr, gdal
except:
    sys.exit('ERROR: cannot find GDAL/OGR modules')

# Set the driver (optional)
wfs_drv = ogr.GetDriverByName('WFS')

# Speeds up querying WFS capabilities for services with alot of layers
gdal.SetConfigOption('OGR_WFS_LOAD_MULTIPLE_LAYER_DEFN', 'NO')

# Set config for paging. Works on WFS 2.0 services and WFS 1.0 and 1.1 with some other services.
gdal.SetConfigOption('OGR_WFS_PAGING_ALLOWED', 'YES')
gdal.SetConfigOption('OGR_WFS_PAGE_SIZE', '10000')

# Open the webservice
url = 'http://example-service.com/wfs'
wfs_ds = wfs_drv.Open('WFS:' + url)
if not wfs_ds:
    sys.exit('ERROR: can not open WFS datasource')
else:
    pass

# iterate over available layers
for i in range(wfs_ds.GetLayerCount()):
    layer = wfs_ds.GetLayerByIndex(i)
    srs = layer.GetSpatialRef()
    print 'Layer: %s, Features: %s, SR: %s...' % (layer.GetName(), layer.GetFeatureCount(), srs.ExportToWkt()[0:50])

    
# set up the shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile")

# create the data source
data_source = driver.CreateDataSource("volcanoes.shp")

# create the spatial reference, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# create the layer
layer = data_source.CreateLayer("volcanoes", srs, ogr.wkbPoint)

# Add the fields we're interested in
field_name = ogr.FieldDefn("Name", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)
field_region = ogr.FieldDefn("Region", ogr.OFTString)
field_region.SetWidth(24)
layer.CreateField(field_region)
layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Elevation", ogr.OFTInteger))

# Process the text file and add the attributes and features to the shapefile
# iterate over features
feat = layer.GetNextFeature()
while feat is not None:
  feat = layer.GetNextFeature()
    
  # create the feature
  feature = ogr.Feature(layer.GetLayerDefn())
  # Set the attributes using the values from the delimited text file
  feature.SetField("Name", feat['Name'])
  feature.SetField("Region", feat['Region'])
  feature.SetField("Latitude", feat['Latitude'])
  feature.SetField("Longitude", feat['Longitude'])
  feature.SetField("Elevation", feat['Elev'])

  # create the WKT for the feature using Python string formatting
  wkt = "POINT(%f %f)" %  (float(feat['Longitude']) , float(feat['Latitude']))

  # Create the point from the Well Known Txt
  point = ogr.CreateGeometryFromWkt(wkt)

  # Set the feature geometry using the point
  feature.SetGeometry(point)
  # Create the feature in the layer (shapefile)
  layer.CreateFeature(feature)
  # Dereference the feature
  feature = None
  feat = None


# Save and close the data source
data_source = None
1

You can do this from the command line using ogr2ogr:

ogr2ogr -of "ESRI Shapefile \
        -spat 419000 419000 421000 42100 \ 
        data_dump.shp
        WFS:https://environment.data.gov.uk/spatialdata/flood-map-for-planning-rivers-and-sea-flood-zone-2/wfs

You can call this command from within Python, or you can use the gdal.VectorTranslate function.

1

The WFS provides the functionality (in theory) to allow you to request the data directly, you just need to write a URL to request the data.

A WFS request would look like:

https://environment.data.gov.uk/spatialdata/flood-map-for-planning-rivers-and-sea-flood-zone-2/wfs?service=WFS&version=2.0.0&request=GetFeature&typeNames=ms:Flood_Map_for_Planning_Rivers_and_Sea_Flood_Zone_2&srsName=urn:ogc:def:crs:EPSG::27700&419000,419000,421000,421000&OUTPUTFORMAT=SHAPE%2BZIP&

in theory because though the service reports it supports <ows:Value>SHAPE+ZIP</ows:Value> as an output format in the GetCapabilities response, the service reports it doesn't support it on doing a GetFeature request; it might be worth contacting the service providers and ask them to fix the issue

You can use owslib to work directly with WFS like:

from owslib.wfs import WebFeatureService

eaWFS = WebFeatureService(url='https://environment.data.gov.uk/spatialdata/flood-map-for-planning-rivers-and-sea-flood-zone-2/wfs?', version='2.0.0')

floodData = eaWFS.getfeature(typename='ms:Flood_Map_for_Planning_Rivers_and_Sea_Flood_Zone_2', bbox=(419000,419000,421000,421000), outputFormat='geopackage')

out = open('C:\\Testing\data.gpkg', 'wb')
out.write(floodData.read())
out.close()

Using geopackage as a substitute for shapefile

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.