1

I have some geoscript code combines WFS and shapefile polygon layers together and writes the result as a GML file. It works fine. However, the client requires a new feature which should combine different geometry types (polygon, polyline and points) together and output them as a GML file. As a geoscript layer's schema limitation, everything put into the layer will be transferred to the schemas' type. For example, if I set the schema as multipolygon and put a point into it, the point will be transferred to a polygon with four same coordinates.

Here is my geoscript code:

from geoserver.wps import process
from geoscript.feature import Feature
from geoscript.layer import Layer
from geoscript.workspace import Directory
from geoscript import geom
from geoscript import proj
from geoscript.feature import Schema
from geoscript.layer.io.gml import readGML
from geoscript.layer.io.gml import writeGML
import os


@process(
        inputs = {
        'bound': (str, 'the boundry to cut shps'),
        'uid': (str, 'the uid of current request')
        },
        outputs = {
        'output1': (str, 'test'),
        'output2': (str, 'test')
        }
)

def run(bound, uid):

        # Set source/destination directories
        wfs_source_dir = '/some/directory/'
        shp_source_dir = Directory('/some/directory/')
        des_dir = '/some/directory/'

        output1 = ''
        output2 = ''

        # Get cutting boundary
        boundary = geom.fromWKT(bound)

        # Get wfs layers and put them in wfs_layers
        wfs_layers = []

        for filename in os.listdir(wfs_source_dir):
                file = open(wfs_source_dir + '/' + filename, 'r')
                layer = readGML(file)
                wfs_layers.append(layer.reproject('epsg:3857'))
        has_wfs = 0
        if len(wfs_layers) > 0:
                has_wfs = 1

        # Get shp layers and put them in shp_layers
        shp_name_list = shp_source_dir.layers()

        shp_layers = []
        for shpName in shp_name_list:
                shp = shp_source_dir.get(shpName)
                shp_layers.append(shp.reproject('epsg:3857'))

        # Create new schema which is combination of all layers
        schema = ''
        fields = []

        if has_wfs == 1:
                for index, layer in enumerate(wfs_layers):
                        if index == 0:
                                fields = layer.schema.fields
                        else:
                                for field in layer.schema.fields[4:]:
                                        fields.append(field)
                for index, layer in enumerate(shp_layers):
                        for field in layer.schema.fields[1:]:
                                fields.append(field)
        else:
                for index, layer in enumerate(shp_layers):
                        if index == 0:
                                fields = layer.schema.fields
                        else:
                                for field in layer.schema.fields[1:]:
                                        fields.append(field)

        schema = Schema(uid.replace('-', ''), fields)
        outlayer.add({'the_geom': geom.Point(-12577404,4914391)})

        # Create new layer and add all features to it
        outlayer = Layer(schema = schema)
        for layer in wfs_layers:
                cursor = layer.cursor()
                features = cursor.read(layer.count())
                for feature in features:
                        geometry = feature.geom
                        if geometry.intersects(boundary):
                                new_feature = feature
                                new_geometry = geometry.intersection(boundary)
                                new_feature.set('the_geom', new_geometry)
                                outlayer.add(new_feature)
                cursor.close()
        for layer in shp_layers:
                cursor = layer.cursor()
                features = cursor.read(layer.count())
                for feature in features:
                        geometry = feature.geom
                        if geometry.intersects(boundary):
                                new_feature = feature
                                new_geometry = geometry.intersection(boundary)
                                new_feature.set('the_geom', new_geometry)
                                outlayer.add(new_feature)
                cursor.close()

        # Write new layer as GML
        out = open(des_dir + '/' + uid + '.gml', 'w')
        writeGML(outlayer, output = out)

        return {'output1': output1, 'output2': output2}
  1. Does GML has capability to handle different geometry types?

  2. Is it possible to add multiple geoscript layers to a single gml file?

  3. If no, what should I do to solve my problem?

  • could you just use a Geometry for the geometry attribute instead of a Polygon? – Ian Turton Jun 29 '17 at 12:17
  • Thank you so much @iant. Your solution solved my issue! I didn't know it accept geom.Geometry. Could you put this in an answer so that I can accept it? – zhangjinzhou Jun 29 '17 at 17:59
1

You should be able to use a type of geometry which is the super type of all the other types.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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