<div dir="ltr"><div>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.</div><div>anyway there is a little issue in the script you provide, even the sample point has to be transformed to MapCanvas crs <br>p1 = transform.trasform(QgsPointXY(1.60301, 38.38453))</div><div><div>In effect putting the project to the same crs of the mesh layer (4326) the script outputs the correct values</div></div><div>so the above script should be the following</div><div><br></div><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)<span class="gmail-im"><br><br>p1 = transform.trasform(QgsPointXY(1.60301, 38.38453))<br><br>ds_group_index = 1<br><br>for g in mesh.datasetGroupsIndexes():<br></span> 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()<span class="gmail-im"><br> value = mesh.datasetValue(dataset, p1)<br> print ("VALUES P1",t,value.scalar(),value.x(),value.y())</span></div><div>##################################################################################</div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno mer 3 feb 2021 alle ore 17:14 Vincent Cloarec <<a href="mailto:vcloarec@gmail.com">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>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>
</blockquote></div>