[QGIS-Developer] Can't get values from QgsMeshLayer

[Vincent Cloarec]

I rewrite your script :

##############################################################
mesh=QgsMeshLayer('/your/path//20210202_211414_ICON_EU_P06_.grb2','mesh','mdal')
print("mesh is valid:",mesh.isValid())

# First, we need to build the triangular mesh used to render and access the
mesh from the layer
# QgsCoordinatesTransform is needed, here we build one with destination to
the map but could be other one
# If the layer is already render in the map, no need to do that if we want
it in the same coordinates than the map

destinationCrs=iface.mapCanvas().mapSettings().destinationCrs()
transform=QgsCoordinateTransform(mesh.crs(),destinationCrs,QgsProject.instance().transformContext())
mesh.updateTriangularMesh(transform)

p1 = QgsPointXY(1.60301, 38.38453)

ds_group_index = 1

for g in mesh.datasetGroupsIndexes():
    print ("MESH GROUP:", g,
mesh.datasetGroupMetadata(QgsMeshDatasetIndex(g,0)).name())

for i in range(mesh.datasetCount(QgsMeshDatasetIndex(ds_group_index,0))):
    dataset = QgsMeshDatasetIndex(ds_group_index, i)
    meta = mesh.datasetMetadata(dataset)
    t = meta.time()
    value = mesh.datasetValue(dataset, p1)
    print ("VALUES P1",t,value.scalar(),value.x(),value.y())
##################################################################################


Be careful to not mix QgsMeshDatasetIndex index from provider with
QgsMeshDatasetIndex from layer. Indeed, since 3.16 we can have extra
dataset group (for example virtual one), the layer handles with its own way
the group indexes

Le mer. 3 févr. 2021 à 11:31, Enrico Ferreguti <enricofer at gmail.com> a
écrit :

> 3.16
> The attachment was missing.
> You can download it from here:
> https://beato.duckdns.org/nextcloud/s/Ra4opS94zYfgfXa
>
> Il mer 3 feb 2021, 16:26 Vincent Cloarec <vcloarec at gmail.com> ha scritto:
>
>> Hi Enrico,
>>
>> As that will depend of the version of QGIS, which version do you use?
>>
>> Le mer. 3 févr. 2021 à 10:59, Enrico Ferreguti <enricofer at gmail.com> a
>> écrit :
>>
>>> I'm trying to write a simple function for retrieving wind speed and
>>> direction at a specified moment from a given grib file readable as mesh
>>> layer the prototype script is the following:
>>>
>>> mesh = QgsMeshLayer("[path to the attached
>>> grib]/20210202_211414_ICON_EU_P06_.grb2","test","mdal")
>>> #mesh = iface.activeLayer()
>>> p1 = QgsPointXY(1.60301, 38.38453)
>>>
>>> ds_group_index = 1
>>>
>>> for g in mesh.datasetGroupsIndexes():
>>>     print ("MESH GROUP:", g,
>>> mesh.datasetGroupMetadata(mesh.datasetIndexAtRelativeTime(h0,g)).name())
>>>
>>> for i in range(mesh.dataProvider().datasetCount(ds_group_index)):
>>>     meta =
>>> mesh.dataProvider().datasetMetadata(QgsMeshDatasetIndex(ds_group_index, i))
>>>     t = meta.time()
>>>     dataset = QgsMeshDatasetIndex(ds_group_index, i)
>>>     value = mesh.datasetValue(dataset, p1)
>>>     print ("VALUES P1",t,value.scalar(),value.x(),value.y())
>>>
>>> When I run the script I can get dataset names from metadata but
>>> I'getting Nan as values
>>> I verified that the grib is loadable in mapcanvas and the sample point
>>> is within grib extent.
>>> I read the note in api (
>>> https://qgis.org/api/classQgsMeshLayer.html#a76a7ec072b4acaf45c1a10f5872c4948)
>>> "It uses previously cached and indexed triangular mesh and so if the layer
>>> has not been rendered previously (e.g. when used in a script) it returns
>>> NaN value" so I tried to process an already loaded mesh (decomment line 2)
>>> without luck.
>>>
>>> Is there anybody who can help me finding what I'm doing wrong? Thanks.
>>> _______________________________________________
>>> QGIS-Developer mailing list
>>> QGIS-Developer at lists.osgeo.org
>>> List info: https://lists.osgeo.org/mailman/listinfo/qgis-developer
>>> Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-developer
>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-developer/attachments/20210203/4fd41c26/attachment.html>