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

[Enrico Ferreguti]

Thanks for the quick reply Vincent, it was really the crs mismatch. I
didn't aware that for extracting values the underlying triangular mesh has
to be updated.
anyway there is a little issue in the script you provide, even the sample
point has to be transformed to MapCanvas crs
p1 = transform.trasform(QgsPointXY(1.60301, 38.38453))
In effect putting the project to the same crs of the mesh layer (4326) the
script outputs the correct values
so the above script should be the following

##############################################################
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 = transform.trasform(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())
##################################################################################

Il giorno mer 3 feb 2021 alle ore 17:14 Vincent Cloarec <vcloarec at gmail.com>
ha scritto:

> 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/7263e49f/attachment.html>