Procesos web geoespaciales:

dificultades y alternativas en el caso del visor SACOSTA

Antes de comenzar

Presentación realizada con reveal.js (HTML presentations made easy)

https://slid.es/bielfrontera/procesos_geoespaciales_sacosta/live

SACOSTA

Visor Sensibilidad Ambiental de la línea de Costa

http://gis.socib.es/sacosta

PostGIS

GeoServer

GeoExplorer

Nueva funcionalidad

Estudio de costa afectada por vertido

Problema geoespacial

A partir de una región dada, ¿de qué manera podemos obtener una estadística del tipo de costa?

Aproximación utilizando

WPS

Documentación GeoServer

How many miles of roads are crossing a protected area?

Geoprocesos de GeoServer /
JTS Topology Suite

gs:IntersectionFeatureCollection

gs:CollectGeometries

JTS:length

XML request


<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">

  <ows:Identifier>JTS:length</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>geom</ows:Identifier>
      <wps:Reference mimeType="text/xml; subtype=gml/3.1.1"
                      xlink:href="http://geoserver/wps" method="POST">
      <wps:Body>
        <wps:Execute version="1.0.0" service="WPS">
          <ows:Identifier>gs:CollectGeometries</ows:Identifier>
          <wps:DataInputs>
            <wps:Input>
              <ows:Identifier>features</ows:Identifier>
              <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wps" method="POST">
                <wps:Body>
                  <wps:Execute version="1.0.0" service="WPS">
                    <ows:Identifier>gs:IntersectionFeatureCollection</ows:Identifier>
                    <wps:DataInputs>
                      <wps:Input>
                        <ows:Identifier>first feature collection</ows:Identifier>
                        <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wfs" method="POST">
                          <wps:Body>
                            <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
                              <wfs:Query typeName="sf:roads"/>
                            </wfs:GetFeature>
                          </wps:Body>
                        </wps:Reference>
                      </wps:Input>
                      <wps:Input>
                        <ows:Identifier>second feature collection</ows:Identifier>
                        <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wfs" method="POST">
                          <wps:Body>
                            <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
                              <wfs:Query typeName="sf:restricted"/>
                            </wfs:GetFeature>
                          </wps:Body>
                        </wps:Reference>
                      </wps:Input>
                      <wps:Input>
                        <ows:Identifier>first attributes to retain</ows:Identifier>
                        <wps:Data>
                          <wps:LiteralData>the_geom cat</wps:LiteralData>
                        </wps:Data>
                      </wps:Input>
                      <wps:Input>
                        <ows:Identifier>second attributes to retain</ows:Identifier>
                        <wps:Data>
                          <wps:LiteralData>cat</wps:LiteralData>
                        </wps:Data>
                      </wps:Input>
                    </wps:DataInputs>
                    <wps:ResponseForm>
                      <wps:RawDataOutput mimeType="text/xml;
                                         subtype=wfs-collection/1.0">
                        <ows:Identifier>result</ows:Identifier>
                      </wps:RawDataOutput>
                    </wps:ResponseForm>
                  </wps:Execute>
                </wps:Body>
              </wps:Reference>
            </wps:Input>
          </wps:DataInputs>
          <wps:ResponseForm>
            <wps:RawDataOutput mimeType="text/xml; subtype=gml/3.1.1">
              <ows:Identifier>result</ows:Identifier>
            </wps:RawDataOutput>
          </wps:ResponseForm>
        </wps:Execute>
      </wps:Body>
    </wps:Reference>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput>
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

Proceso 1

 <wps:Execute version="1.0.0" service="WPS">
  <ows:Identifier>gs:IntersectionFeatureCollection</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
        <wfs:Query typeName="sf:roads"/>
      </wfs:GetFeature>
    </wps:Input>
    <wps:Input>
      <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
        <wfs:Query typeName="sf:restricted"/>
      </wfs:GetFeature>
    </wps:Input>
  </wps:DataInputs>
</wps:Execute>

Proceso 2

 <wps:Execute version="1.0.0" service="WPS">
  <ows:Identifier>gs:CollectGeometries</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>features</ows:Identifier>
      <wps:Reference mimeType="text/xml; subtype=wfs-collection/1.0" xlink:href="http://geoserver/wps" method="POST">
        <wps:Body>
          <wps:Execute>
            ________________ PROCESO 1 ________________
          </wps:Execute> 
        </wps:Body>
      </wps:Reference>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="text/xml; subtype=gml/3.1.1">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

Proceso 3

 <wps:Execute version="1.0.0">
  <ows:Identifier>JTS:length</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>geom</ows:Identifier>
      <wps:Reference mimeType="text/xml; subtype=gml/3.1.1"
                      xlink:href="http://geoserver/wps" method="POST">
      <wps:Body>
        <wps:Execute version="1.0.0" service="WPS">
          ________________ PROCESO 2 ________________
        </wps:Execute>
      </wps:Body>
    </wps:Reference>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput>
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

En nuestro caso

No encontramos una manera sencilla de devolver longitud por tipo utilizando los procesos ya definidos
Tendríamos que implementar un proceso ad hoc. Reinstalación de GeoServer desarrollo de nuevo proceso (definido como extensión de org.geotools.process.gs.GSProcess )
Podríamos utilizar el cliente de WPS de OpenLayers

Solución ad hoc:

PostGIS + Python Flask

Funciones geoespaciales

en PostGIS

ST_GeomFromText
ST_Transform
ST_Intersects | ST_Intersection
ST_Lenght

Instrucción SQL

SELECT sci."ESICOSTES",
       count(1) as num_features,
       sum(st_length(the_geom_intersec)) as longitud_intersec,
       string_agg(sci."HOTLINK", '|') as hotlink
FROM (
    SELECT sc.*,
    ST_Intersection(
        ST_Transform(ST_GeomFromText('POLYGON ((2.5488447287659 39.528537630451, 2.5492738822082 39.522678517902, 2.5845932105163 39.528636933183, 2.5818037131409 39.534396248632, 2.5488447287659 39.528537630451))', 4326), 3043),
        the_geom) as the_geom_intersec
    FROM sacosta.bal_sa_costa_2012 as sc
    WHERE ST_Intersects(
        ST_Transform(ST_GeomFromText('POLYGON ((2.5488447287659 39.528537630451, 2.5492738822082 39.522678517902, 2.5845932105163 39.528636933183, 2.5818037131409 39.534396248632, 2.5488447287659 39.528537630451))', 4326), 3043),
        the_geom)
) AS sci
GROUP BY "ESICOSTES"
ORDER BY "ESICOSTES"

ST_GeomFromText

ST_GeomFromText('POLYGON ((2.5488447287659 39.528537630451, 2.5492738822082 39.522678517902, 2.5845932105163 39.528636933183, 2.5818037131409 39.534396248632, 2.5488447287659 39.528537630451))', 4326) as selected_region

ST_Transform

ST_Transform(selected_region, 3043) as selected_region_transformed

ST_Intersection

ST_Intersection(  selected_region_transformed,   the_geom) as the_geom_intersec

ST_Length

sum(st_length(the_geom_intersec))

Voilà

Flask

mini framework web de python

@app.route

@app.route('/api/v1.0/sacosta/<polygon_text>', methods=['GET'])
@app.route('/api/v1.0/sacosta/', methods=['POST'])
@crossdomain(origin='*')
def get_sacosta_polygon(polygon_text=None):
    if request.method == 'POST':
        polygon_text = request.form['polygon']

    try:
        polygon = utils.create_polygon(polygon_text)
    except ValueError, e:
        return jsonify({'error': str(e)})
    return utils.send_json_data(gisdata.get_data_sacosta(current_app.config, polygon))

Shapely (Polygon y wkt)

 def create_polygon(polygon_text):
    """ Get a shapely Polygon from text

    :param polygon: string with a list of floats separated by commas that represents
                    polygon vertexs or a WKT representing a polygon
    :returns: Polygon
    """
    if polygon_text.startswith('POLYGON'):
        try:
            # Check polygon is a valid wkt
            polygon = wkt.loads(polygon_text)
        except:
            raise ValueError('Could not create geometry')
    else:
        # get polygon WKT from a list of vertices
        vertices = get_polygon_vertices_from_text(polygon_text)
        polygon = Polygon(vertices)

    return polygon

Psycopg2


def get_data_sacosta(config, polygon):
    """ Get aggregated data of coastline sensibility for a given region

    :param config: app config
    :param polygon: shapely polygon
    :returns: array with the result of longitude and images for each type of coast
    """
    region = "ST_GeomFromText('{polygon_text}', 4326)".format(polygon_text=polygon.wkt)
    sql_sacosta = """... SQL {region} ...""".format(region=region)
    cur = conn.cursor(cursor_factory=DictCursor)
    cur.execute(sql_sacosta)
    # Process the result
    result = []
    for row in cur:
        obj = {
            'esicostes': row['ESICOSTES'],
            'longitud': round(row['longitud_intersec'], 2),
            'hotlink': []
        }
        if row['hotlink'] is not None:
            obj['hotlink'] = [link for link in row['hotlink'].split('|')]
        result.append(obj)
    return result

Despliegue

git
virtualenv
apache2 (pero podría ser nginx)
supervisor
gunicorn (pero podría ser uwsgi)

Código disponible en:

https://github.com/socib/sacosta_gisservices

Integración en GeoExplorer

Nueva herramienta, gxp.plugins.Tool


Ext.namespace("gxp.plugins");
gxp.plugins.SensibilidadAmbiental = Ext.extend(gxp.plugins.Tool, {
    ptype: "sacosta_sensibilidadambiental",
    menuText: "Sensibilidad Ambiental",
    toolTip: "Muestra la sensibilidad ambiental de una área seleccionada",
    services_host: 'gisservices.socib.es',
    constructor: function(config) {
        gxp.plugins.SensibilidadAmbiental.superclass.constructor.apply(this, arguments);
    },
    addActions: function() {
        var actions = gxp.plugins.SensibilidadAmbiental.superclass.addActions.apply(this, [{
            iconCls: "sacosta-icon-sensibilidadambiental",
            enableToggle: true,
            handler: function() {
                this.drawPolygon();
            },
            scope: this,
            toggleHandler: function(something, state) {
            }
        }]);
        return actions;
    }
}

Añadir a la barra de herramientas


GeoExplorer.Composer = {
    config.tools = [
        {...},
        {
            ptype: "sacosta_sensibilidadambiental",
            actionTarget: {
                target: "paneltbar",
                index: 18
            }
        }
    ]
}

Controles OpenLayers

OpenLayers.Control.DrawFeature

 
this.mapControls['draw'] = new OpenLayers.Control.DrawFeature(this.polygonLayer, OpenLayers.Handler.Polygon, {
    'displayClass': 'olControlDrawFeaturePolygon',
    'featureAdded': OpenLayers.Function.bind(this.onFinishDrawPolygon, this)
});

OpenLayers.Control.DragFeature

 
this.mapControls['drag'] = new OpenLayers.Control.DragFeature(this.polygonLayer, {
    'onComplete': OpenLayers.Function.bind(this.getSensibilidadAmbientalInfo, this)
});

Llamada a servicio datos


    getPolygonWKT: function(){
        var wkt_formater = new OpenLayers.Format.WKT({
            'internalProjection': this.target.mapPanel.map.getProjectionObject(),
            'externalProjection': new OpenLayers.Projection("EPSG:4326")
        });
        return wkt_formater.write(this.polygonLayer.features[0]);
    },
    getSensibilidadAmbientalInfo: function() {
        var polygon_wkt = this.getPolygonWKT();

        $.getJSON('http://' + this.services_host + '/api/v1.0/sacosta/' + polygon_wkt, (function(jsondata) {
            if (jsondata.data) {
                this.plotSensibilidadAmbientalInfo(jsondata.data);
            } else {
                console.log('Error retreiving data from gisservices ' + jsondata.error || '');
            }
        }).bind(this));

    },

Gráfico de barras con D3.js

plotSensibilidadAmbientalInfo: function(data) {
    this.svg.selectAll("rect")
        .data(data)
        .enter().append("rect")
        .attr("fill", this.getColorBar.bind(this))
        .attr("x", function(d) {
            return x(d.esicostes);
        })
        .attr("width", x.rangeBand())
        .attr("y", function(d) {
            return y(d.longitud);
        })
        .attr("height", function(d) {
            return height - y(d.longitud);
        });
}

Filtro por tipo de costa


onMouseOverBar: function(d, i) {
    ol_filter = new OpenLayers.Filter.Logical({
        type: OpenLayers.Filter.Logical.AND,
        filters: [
            new OpenLayers.Filter.Comparison({
                type: OpenLayers.Filter.Comparison.EQUAL_TO,
                property: "ESICOSTES",
                value: d.esicostes
            })
        ]
    });
    var filter_1_0 = new OpenLayers.Format.Filter({
        version: "1.0.0"
    });
    var xml = new OpenLayers.Format.XML();
    var new_filter = xml.write(filter_1_0.write(ol_filter));
    layer.params['FILTER'] = new_filter;
    layer.redraw();
};