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

[Vincent Cloarec]

thinking about, if you use only:
mesh.updateTriangularMesh()

no need to handle with crs ans transform, I think results will be in layer
coordinates.


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

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