Moviendo bounding box de imagen, basada en el centroide, hasta la posición contigua con PyQGIS

En post anteriores (1, 2) se consideró la conexión al servicio de Google Satellite y la automatización en un plugin con el fin de obtener tiles de alta resolución georreferenciados. En este post se tratará uno de los algoritmos para que ello fuese posible.

Considera como plantilla la primera imagen de alta resolución que fue grabada a escala 1:2000. A continuación se puede observar el ráster una vez cargado que, por defecto, tiene un Zoom a la capa y que por las características particulares del Map Canvas de mi QGIS lo refleja con una escala de 1:2100.

No obstante, si fijamos manualmente la escala a 1:2000 se puede observar en la imagen siguiente que se ajusta perfectamente al área exhibida por el Map Canvas.

Sin embargo, no importa el nivel inicial de Zoom de la capa. Lo necesario es que ésta sea simplemente capa activa. El código a continuación lo primero que hace es crear un objeto de la clase QgsRectangle para determinar el ancho (que se usará en la translación) y la geometría (a partir del formato WKT) de la bounding box. Esta última servirá para determinar el centroide y las nuevas coordenadas resultantes del desplazamiento. La función afín dará cuenta de la translación de los puntos de la bounding box que servirán de base para crear una memory layer tipo polígono y empalmar de la manera adecuada.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

layer = iface.activeLayer() extent = layer.extent() width = extent.width() pol_extent = extent.asWktPolygon() print pol_extent geom_pol = QgsGeometry.fromWkt(pol_extent) geom = geom_pol.asPolygon() geom_t = geom centroid = geom_pol.centroid() c1 = centroid.asPoint() c2 = QgsPoint(c1.x() + width, c1.y()) #calculating displacement based in direction vector px = c2.x() - c1.x() py = c2.y() - c1.y() #translation as an affine transformation for i in range(len(geom_t)): n = len(geom_t[i]) for j in range(n): geom_t[i][j].setX(geom[i][j].x() + px) geom_t[i][j].setY(geom[i][j].y() + py) d = QgsGeometry.fromPolygon(geom_t).exportToWkt() print d epsg = layer.crs().postgisSrid() uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes" mem_layer = QgsVectorLayer(uri, 'bounding_box', 'memory') prov = mem_layer.dataProvider() feat = QgsFeature() feat.setAttributes([0]) feat.setGeometry(QgsGeometry.fromWkt(d)) prov.addFeatures([feat]) QgsMapLayerRegistry.instance().addMapLayer(mem_layer) iface.zoomToActiveLayer() mapcanvas = iface.mapCanvas() mapcanvas.zoomScale(2000)

Una vez ejecutado el código anterior se obtiene el resultado siguiente:

Observe que la Python Console no está empotrada en el Map Canvas para aprovechar su superficie máxima y que ésta sea equivalente a la de la imagen plantilla cuando se aplique la escala 1:2000. Al conectar el servicio de Google Satellite se puede observar en la imagen a continuación el tile correspondiente:

Para corroborar que es el tile esperado, como en la Python Console se imprime también la geometría en formato WKT para el primer polígono, las colocamos juntas en la imagen siguiente (Zoom Full) con un 50 % de transparencia. Efectivamente son tiles contiguos.

Nota:

Este códido funciona en QGIS3.

1 2 3 4

import requests service_url = "mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}" service_uri = "type=xyz&zmin=0&zmax=21&url=https://"+requests.utils.quote(service_url) tms_layer = iface.addRasterLayer(service_uri, "Google Sat", "wms")