WFS layer saved as a shapefile

Question

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

Answer 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

Answer 2

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.

Answer 3

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