Using the Intersect filter in OWSlib to filter data between two WFS

Question

I'm trying to use OWSLib to deal with multiple WSF features in python, for now, I'm trying to return all features in sfb WFS that intersects with the specifc feature retrieved in sicar WFS.

from owslib.wfs import WebFeatureService, Authentication
from owslib.fes2 import PropertyIsEqualTo, Intersects
from owslib.etree import etree

auth = Authentication()
auth.verify = False

sicar = WebFeatureService(
    url="https://geoserver.car.gov.br/geoserver/sicar/ows", auth=auth, version='2.0.0')

sfb = WebFeatureService(
    'http://sistemas.florestal.gov.br/geoserver/ows', version='2.0.0')

car = 'MT-5103700-92A39DC62FDF4D88874F218D41BE1C7F'

filter_car = PropertyIsEqualTo('cod_imovel', car).toXML()

carInfo = sicar.getfeature(typename='sicar:sicar_imoveis_mt',
                           filter=filter_car, method='post', outputFormat='json')

The code above is working fine, the problem starts when i try to get the Intersects

geom_property = 'geo_area_imovel';
carInfoStr = carInfo.read().decode('utf-8');
geom = json.loads(carInfoStr)['features'][0]['geometry']['coordinates']

#I'm trying to pass the the coordinates list to Intersects, but i tried with other keys of the object and with the xml response, but doesn't wotk as well.

filter_contains = Intersects(geom_property, geom).toXML()

intersected = sfb.getfeature(typename='CNFP_orig:APP_Total_simp',
                       filter=filter_contains, outputFormat='json', method='post')

The returned errors involve the toXML function, wich depends on the geometry i pass to the intersects, in this case its a string, but need to be the constructor of something of the owslib, but the lib only have Point constructor.

I could pass a point if I could retrieve the centroid of geometry, but I don't now how to do this with this lib. BUT, even the example passed in a official PR is not working.

the PR with example: https://github.com/geopython/OWSLib/pull/780

The error in this code:

Traceback (most recent call last):
  File "/home/gabrielubuntu/repos/pygeogis/app.py", line 38, in <module>
    filter_contains = Intersects(geom_property, geom).toXML()
  File "/home/gabrielubuntu/repos/pygeogis/.venv/lib/python3.10/site-packages/owslib/fes2.py", line 420, in toXML
    node.append(self.geometry.toXML())
AttributeError: 'dict' object has no attribute 'toXML'

Answer 1

Strangely this seems to be very hard, so I've probably missed something, but the following code seems to work:

import lxml.etree as ET                                                                                                 
from owslib.fes2 import PropertyIsEqualTo                                                                               
from owslib.wfs import WebFeatureService, Authentication                                                                
                                                                                                                        
auth = Authentication()                                                                                                 
auth.verify = False                                                                                                     
                                                                                                                        
sicar = WebFeatureService(                                                                                              
    url="https://geoserver.car.gov.br/geoserver/sicar/ows",                                                             
    auth=auth, version='2.0.0')                                                                                         
                                                                                                                        
                                                                                                                        
car = 'MT-5103700-92A39DC62FDF4D88874F218D41BE1C7F'                                                                     
# Here the filter works as expected                                                                                                                        
filter_car = PropertyIsEqualTo('cod_imovel', literal=car)                                                               
                                                                                                                      
print(filter_car)                                                                                                       
print("--"*10)                                                                                                          
                                                                                                                        
carInfo = sicar.getfeature(typename='sicar:sicar_imoveis_mt',                                                           
                           filter=filter_car.toXML(), method='post')                                                    

# should probably look this up as below                                                                                                                    
geom_property = 'geo_area_imovel'      
# if you prefer JSON I'm sure that would work too                                                                                 
cars = ET.parse(carInfo)                                                                                                
ET.dump(cars.getroot())                                                                                                 
print("--"*10)                                                                                                          
                                                                                                                        
geom = cars.find(".//{geoserver.car.gov.br}"+geom_property)                                                             
sfb = WebFeatureService('http://sistemas.florestal.gov.br/geoserver/ows',                                               
                        version='2.0.0')                                                                                
name = 'CNFP_orig:APP_Total_simp'                                                                                       
schema = sfb.get_schema(name)                                                                                           
print(schema)                                                                                                           
print("--"*10)  

# I can't see how you can use the Intersects object here as I can't 
# work out how to get a geometry object from the GML, so I build an XML
# Filter (we could add namespaces but GeoServer is lazy about that)                                                                                                        
filter = ET.Element("Filter")                                                                                           
intersects = ET.SubElement(filter, 'Intersects')                                                                        
propName = ET.SubElement(intersects, 'PropertyName')                                                                    
propName.text = schema.get('geometry_column')                                                                           
intersects.append(geom)                                                                                                 
                                                                                                                        
ET.dump(filter)                                                                                                         
print("--"*10)                                                                                                          

# I couldn't make this work with a POST - so this might fail with very
# large geometries                                                                                                                        
intersected = sfb.getfeature(typename=name,                                                                             
                             filter=filter)                                                                             
                                                                                                                        
result = ET.parse(intersected)                                                                                          
ET.dump(result.getroot())