3

I'm trying to open a raster image of 15x15Km with a google basemap with PyQGIS. When the code finishes, I want the map canvas to zoom into the raster image. It almost works and I can add rasters and even vectors, but if I add a basemap, it's like the basemap overides the zoom level. so this simple code:

#basemap
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")

#image
image=”/media/patrice/DataDrive/GORiM/ManualFineTuning/GEE/S2_S025_Z12R.tif”
Cname=os.path.basename(image)[:-4]
Zone=os.path.basename(image)[-8:-4]
Sample=os.path.basename(image)[3:7]
Cname=os.path.basename(image)[:-4]
rlayer = iface.addRasterLayer(image,Cname,"gdal")
vLayer = iface.activeLayer()
extent = vLayer.extent()
iface.mapCanvas().setExtent(extent)
iface.mapCanvas().refresh()

Loads the layer and the basemap but the display has the whole world, in the correct projection, but fully zoomed out.

I also tried to get the active layer by name:

Cname=os.path.basename(image)[:-4]
rlayer = iface.addRasterLayer(image,Cname,"gdal")
layer =  QgsProject.instance().mapLayersByName(Cname)[0]
layer.loadNamedStyle('/media/patrice/DataDrive/GORiM/ManualFineTuning/ImageNIRStyle.qml')
layer.triggerRepaint()
layer =  QgsProject.instance().mapLayersByName(Cname)[0]
extent = layer.extent()
iface.mapCanvas().setExtent(extent)
iface.mapCanvas().refresh()

With the same result and no error messages. I think the problem is in setting the active layer, but I have tried printing layer names and they match.

UPDATE
No luck with the answer below. Both snipets of code gave the same result. I also investigated the idea that it was a conflict between layer CRS and project crs by adding:

Cname = Path(image).stem
Zone = Cname.split('_')[-1]
crs=str(GetCRS(Zone[-1], Zone[-3:-1]))
QGIScrs=QgsCoordinateReferenceSystem(crs)
QgsProject.instance().setCrs(QGIScrs) 

where GetCRS gets the EPSG code from the zone number. The zones seem correct and the project crs matches the layer crs. Also, the google basemap is centered on the correct spot and clearly not in 3857 or in 4326. So is this an actual bug?

Raster File Raster file sample below. It's a 3 band Sentinel 2 file generated by Google Earth Engine. If this matters, I'm using Ubuntu 20.04 and 3.22.

https://drive.google.com/file/d/1wakV7g_eYNjq0NFv27zAq_FUVvq7heZe/view?usp=sharing

3
  • Why not use the extent of your rlayer? If I understand you correctly, that is what you want. :D Commented Dec 9, 2022 at 13:22
  • Tried that but same results. I've also printed the rlayer (and others) and it prints the correct layer name. Just won't change them. Commented Dec 9, 2022 at 13:38
  • Please see update at the end of my answer.
    – Ben W
    Commented Dec 14, 2022 at 2:40

1 Answer 1

6

I suspect that the problem is a mis-match between your project CRS and your raster CRS. Bear in mind that the extent() method of QgsMapLayer returns a QgsRectangle in the layer CRS. When you pass a QgsRectangle object to QgsMapCanvas.setExtent(), it must be in the same CRS as the canvas.

To work around this in a simple way, I would recommend using iface.zoomToActiveLayer() instead.

The following script should work as expected:

from pathlib import Path
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")

image="/media/patrice/DataDrive/GORiM/ManualFineTuning/GEE/S2_S025_Z12R.tif"

Cname = Path(image).stem
Zone = Cname.split('_')[-1]

rlayer = iface.addRasterLayer(image, Cname, "gdal")
rlayer.loadNamedStyle("/media/patrice/DataDrive/GORiM/ManualFineTuning/ImageNIRStyle.qml")
rlayer.triggerRepaint()

iface.setActiveLayer(rlayer)
iface.zoomToActiveLayer()

Otherwise, if you want to use iface.mapCanvas().setExtent() you will need to transform the extent from the layer CRS to the project CRS. To do that you can use the QgsCoordinateTransform class.

E.g.

from pathlib import Path
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")

image="/media/patrice/DataDrive/GORiM/ManualFineTuning/GEE/S2_S025_Z12R.tif"

Cname = Path(image).stem
Zone = Cname.split('_')[-1]

rlayer = iface.addRasterLayer(image, Cname, "gdal")
rlayer.loadNamedStyle("/media/patrice/DataDrive/GORiM/ManualFineTuning/ImageNIRStyle.qml")
rlayer.triggerRepaint()

def transform(rect):
    rlayer_crs = rlayer.crs()
    canvas_crs = QgsProject.instance().crs()
    if rlayer_crs != canvas_crs:
        xform = QgsCoordinateTransform(rlayer_crs, canvas_crs, QgsProject.instance())
        rect = xform.transform(rect)
        enlargement = min(rect.width(), rect.height())*0.05
    return rect.buffered(enlargement)
    
iface.mapCanvas().setExtent(transform(rlayer.extent()))

EDIT- After further testing: While both code snippets above worked in 3.28 on Ubuntu 22.04, I found a similar behaviour to that described by the OP in 3.22 on Windows though only when running the code in a new, empty project.

A workaround is to load the layers in two separate functions and call the second function from the first, using QTimer.singleShot to apply a small delay.

from pathlib import Path
from PyQt5.QtCore import QTimer
import requests


def load_raster():
    image="/media/patrice/DataDrive/GORiM/ManualFineTuning/GEE/S2_S025_Z12R.tif"
    Cname = Path(image).stem
    Zone = Cname.split('_')[-1]
    rlayer = iface.addRasterLayer(image, Cname, "gdal")
    iface.setActiveLayer(rlayer)
    iface.zoomToActiveLayer()
    
def load_basemap():
    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")
    if len(QgsProject.instance().mapLayers())<2:
        QTimer.singleShot(100, load_raster)
    else:
        load_raster()

load_basemap()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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