More

Point attribute from underlying raster data TIFFTAG

Point attribute from underlying raster data TIFFTAG


I have a layer of manually set points in qgis that have the following attributes at the moment :

id, lat, long

Additionally, I have a number of different overlapping raster images (sattelite data) that are all tagged with a date using the TIFFTAG_DATETIME in the raster's metadata.

What I would like to do is to form an additional attribute for the point layer that is called dates:

dates = ["2015-01-01", "2015-12-01"]

The list of dates should contain all the dates of satellite imagery that contain the particular point.

Any ideas how to do this in qgis. If needed I can also use different methods of tagging the GeoTiff file. Or maybe even create polygons that have the date as attributes. But the code that is creating and tagging the GeoTiff files runs outside of qgis.


One possibility is to bring the rasters into Postgresql (with postgis extension) and simply selecting the names of the rasters which overlap the point. For python you could install thepsycopg2libraries to do this. To do the import I use theraster2pgsqltool which generates sql which can be fed intopsql.

For more info, see the PostGIS raster page.

Postgis raster offers the following raster functions for bounding boxes and convex hulls…

ST_Envelope() ST_ConvexHull() ST_MinConvexHull()

The last one -ST_MinConvexHull- finds the convex hull without NODATA pixels, leaving you with the polygon outline of data pixels. This makes it suitable for irregularly shaped swathes or surveys.

There's alsoST_Polygon, which returns a multipolygon of data coverage, so I imagine this will allow holes/voids, but i've not been able to check this. Might be worth investigating for DEMs, though.

The possible problem with using QgsRasterLayer.extent() is that returns a bounding box, not a convex hull (or envelope) around the data pixels.

So if you use extent(), your test also needs to also include a test on data/nodata, as Detlev pointed out. I've not tried it with rotated rectangular rasters - but it does work on rectangular rasters with a rotated/irregular data area, surrounded by NODATA.

So assuming you've loaded the rasters into postgres (and adding a timestamp field to each raster) you could do something like this for each point

select timestamp from rastertable where st_contains( st_minconvexhull(rast), st_setsrid( st_makepoint($x,$y), $SRID ) );

Replace$x,$ywith your point coordinate and$SRIDwith the raster SRID

I'm not sure if POSTGIS can extract the TIFFTAG_DATETIME for you, I think this is discarded on import (see this answer) so you need to add this yourself when importing.

But a call togdalinfocertainly can get the tag if you're automating this…

gdalinfo filename.tiff | grep TIFFTAG_DATETIME

You have a point layer, which is the active layer, and a bunch of raster layers from various satellite images, all tagged with TIFFTAG_DATETIME.

I suppose that the images have the same coordinate system and are rectified, no need to determine whether the points are outside the image info. If this should be considered, the script must be extended to identify the image value at the point's position and to decide if the value is data or nodata.

from PyQt4.QtCore import QVariant from osgeo import gdal # dictionary to hold metadata from the images metadata = {} # get point layer and add attribute 'dates' if it does not exist points = iface.activeLayer() prov = points.dataProvider() fni = prov.fieldNameIndex('dates') if fni == -1: if prov.capabilities() & QgsVectorDataProvider.AddAttributes: prov.addAttributes([QgsField('dates', QVariant.String)]) # populate a dict with date tag and bounding box of all raster layers for image_name, image in QgsMapLayerRegistry.instance().mapLayers().items(): if isinstance(image, QgsRasterLayer): # open image with gdal to read metadata ras = gdal.Open(image.dataProvider().dataSourceUri()) datetag = ras.GetMetadata()['TIFFTAG_SOFTWARE'] box = QgsGeometry.fromWkt(layer.extent().asWktPolygon()) metadata[image.name()] = {'datetag': datetag, 'box': box} ras = None # switch to edit mode points.startEditing() # for each point check if it is in any image for feat in points.getFeatures(): dates = [] for image in metadata: result = metadata[image]['box'].contains(feat.geometry()) # if image contains the point add the date to the list if result: dates.append(metadata[image]['datetag']) # cast the list to str and remove brackets feat['dates'] = str(dates).strip('[]') points.updateFeature(feat) # finally save changes points.commitChanges()

I dont think that your problem can be solved with the QGis Gui Interface. I suggest a script like this. It uses the Python Image Processing library to read the desired TiffTag. It is required that all your raster files are loaded in the TOC. It is also required, that you have added a new attribute to your point layer which is named with 'dtattr', type 'string' and maximum length.

##pointsUri=vector from PIL import Image from qgis.utils import * from qgis.core import * points = processing.getObjectFromUri(pointsUri) #get all raster layers form the legendInterface rasterlayers = [layer for layer in iface.legendInterface().layers() if layer.type() == 1] for feature in points.getFeatures(): #check if raster-extend-geometry contains feature geometry datevalues = [] #empty list for the datevalues to store at for rasterlayer in rasterlayers: #create the geometry of the raster extent rasterRectGeom = QgsGeometry.fromWkt(rasterlayer.extent().asWktPolygon()) if rasterRectGeom.contains(feature.geometry()): satImage = Image.open(rasterlayer.source()) tifftags = satImage.tag datevalues.append(tifftags.get(306)) #306 for the TIFFTAG_DATETIME #write to the attribute field points.startEditing() datetimestr = ';'.join(datevalues) #datetimestr =".join(e for e in datevalues) feature['dtattr'] = datetimestr points.updateFeature(feature) points.commitChanges()

You may have to adjust the output a Little bit, as it comes in the yyyy:mm:dd HH:MM:SS Format. I tested the script under QGis 2.8.3 with some example data. To execute the script, create a new 'user script' in the processing toolbar.


Watch the video: arcgis extract selected features. arcgis extract raster values