[QGIS-it-user] Fwd: intersezione layer e calcolo area e crs differenti

[SC]

Salve ho scritto un piccolo script per la console di QGIS per calcolare

le aree intersecate dei vari geometrie di un livello rispetto l'altro, ai piu non sfuggirà la possibilità di usarlo per per calcolare  certificati di destinazione urbanistica rispetto un prg

c'è un problema quando i layer sono  due sistemi di riferimenti
differenti ,il calcolo non funziona, mentre va bene se i liveli sono
sullo stesso crs

sapete dirmi se è

1)un bug di qgis

2) scelta  progettuale (perchè?) da ovviare tramite codice

3)api sbagliate che ho usato si poteva scegliere altra soluzione

saluti



# /* copyright 2017 by Salvo Caligiore caligiore a gmail.com

livelli=iface.mapCanvas().layers()


def dataGeom(geom,lv):
#calcola le aree delle intersezioni delle features di un determinato
livello rispetto ad una data geometria e restituisce una lista di
feature con il corrispondente valore di area intersecato
#la feauture in lista serve per altri calcoli nella procedura principale
                 rect = geom.boundingBox()

                 request=QgsFeatureRequest(rect)
                 lista=[]
                 feat=lv.getFeatures(request)
                 da = QgsDistanceArea()
                 for f in feat:

                     intersezione=f.geometry().intersection(geom)

                     if intersezione==None : continue
                     else :
a=round(da.convertAreaMeasurement(da.measureArea(intersezione),1),2)



                         if a>0 : lista.append((f,round(a,2)))
                         if lv.geometryType()==0: lista.append((f," "))
                 return lista


da = QgsDistanceArea()
#incrocia i livelli fra di loro e interseca ogni feature con le altre
del livello incrociato stampando le informazioni
for x in livelli :
     for y in livelli:
         for f in y.getFeatures():

             lst=dataGeom(f.geometry(),x)

             for z,a in  lst:
areafeature=round(da.convertAreaMeasurement(f.geometry().area(),1),2)
                 print x.name(), y.name(),"ID",f.id(),"Area
Totale",areafeature , "Area Intersecata", a , "ZID",z.id(),"rap % fra
inter e area f" ,str(round((a/areafeature)*100,2))+'%'