<div dir="ltr"><div>I rewrite your script :</div><div><br></div><div>##############################################################<br></div><div>mesh=QgsMeshLayer('/your/path//20210202_211414_ICON_EU_P06_.grb2','mesh','mdal')<br>print("mesh is valid:",mesh.isValid())<br><br># First, we need to build the triangular mesh used to render and access the mesh from the layer<br># QgsCoordinatesTransform is needed, here we build one with destination to the map but could be other one</div><div># 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</div><div><br></div><div>destinationCrs=iface.mapCanvas().mapSettings().destinationCrs()<br>transform=QgsCoordinateTransform(mesh.crs(),destinationCrs,QgsProject.instance().transformContext())<br>mesh.updateTriangularMesh(transform)<br><br>p1 = QgsPointXY(1.60301, 38.38453)<br><br>ds_group_index = 1<br><br>for g in mesh.datasetGroupsIndexes():<br> print ("MESH GROUP:", g, mesh.datasetGroupMetadata(QgsMeshDatasetIndex(g,0)).name())<br><br>for i in range(mesh.datasetCount(QgsMeshDatasetIndex(ds_group_index,0))):<br> dataset = QgsMeshDatasetIndex(ds_group_index, i)<br> meta = mesh.datasetMetadata(dataset)<br> t = meta.time()<br> value = mesh.datasetValue(dataset, p1)<br> print ("VALUES P1",t,value.scalar(),value.x(),value.y())</div><div>##################################################################################</div><div><br></div><div><br></div><div>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<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le mer. 3 févr. 2021 à 11:31, Enrico Ferreguti <<a href="mailto:enricofer@gmail.com" target="_blank">enricofer@gmail.com</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto">3.16<div dir="auto">The attachment was missing. </div><div dir="auto">You can download it from here: <a href="https://beato.duckdns.org/nextcloud/s/Ra4opS94zYfgfXa" target="_blank">https://beato.duckdns.org/nextcloud/s/Ra4opS94zYfgfXa</a></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il mer 3 feb 2021, 16:26 Vincent Cloarec <<a href="mailto:vcloarec@gmail.com" target="_blank">vcloarec@gmail.com</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi Enrico,</div><div><br></div><div>As that will depend of the version of QGIS, which version do you use?<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le mer. 3 févr. 2021 à 10:59, Enrico Ferreguti <<a href="mailto:enricofer@gmail.com" rel="noreferrer" target="_blank">enricofer@gmail.com</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>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:</div><div><br></div><div>mesh = QgsMeshLayer("[path to the attached grib]/20210202_211414_ICON_EU_P06_.grb2","test","mdal")</div><div>#mesh = iface.activeLayer()<br></div><div>p1 = QgsPointXY(1.60301, 38.38453)<br><br>ds_group_index = 1<br><br>for g in mesh.datasetGroupsIndexes():<br> print ("MESH GROUP:", g, mesh.datasetGroupMetadata(mesh.datasetIndexAtRelativeTime(h0,g)).name())<br><br>for i in range(mesh.dataProvider().datasetCount(ds_group_index)):<br> meta = mesh.dataProvider().datasetMetadata(QgsMeshDatasetIndex(ds_group_index, i))<br> t = meta.time()<br> dataset = QgsMeshDatasetIndex(ds_group_index, i)<br> value = mesh.datasetValue(dataset, p1)<br> print ("VALUES P1",t,value.scalar(),value.x(),value.y())</div><div><br></div><div>When I run the script I can get dataset names from metadata but I'getting Nan as values</div><div>I verified that the grib is loadable in mapcanvas and the sample point is within grib extent.</div><div>I read the note in api (<a href="https://qgis.org/api/classQgsMeshLayer.html#a76a7ec072b4acaf45c1a10f5872c4948" rel="noreferrer" target="_blank">https://qgis.org/api/classQgsMeshLayer.html#a76a7ec072b4acaf45c1a10f5872c4948</a>) "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.</div><div><br></div><div>Is there anybody who can help me finding what I'm doing wrong? Thanks.<br></div></div>
_______________________________________________<br>
QGIS-Developer mailing list<br>
<a href="mailto:QGIS-Developer@lists.osgeo.org" rel="noreferrer" target="_blank">QGIS-Developer@lists.osgeo.org</a><br>
List info: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-developer" rel="noreferrer noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-developer</a><br>
Unsubscribe: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-developer" rel="noreferrer noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/qgis-developer</a><br>
</blockquote></div>
</blockquote></div>
</blockquote></div>