4

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.

2
  • AFAIK the point sampling tool has the nearest functionality what you are looking for, but can't solve your problem. Is it really good to put several dates into one field in the database? I suppose you had better to write a Python script using GDAL to solve your task.
    – Zoltan
    Commented Dec 7, 2015 at 8:53
  • Thanks. My thought so far was to create a polygon for each satellite image and tag it with it's creation date. Then I could use a method like. "Join attributes by location" to extract the dates from the polygon layer to the point layer. Commented Dec 7, 2015 at 12:30

3 Answers 3

3
+50

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 the psycopg2 libraries to do this. To do the import I use the raster2pgsql tool which generates sql which can be fed into psql.

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 also ST_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, $y with your point coordinate and $SRID with 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 to gdalinfo certainly can get the tag if you're automating this...

gdalinfo filename.tiff | grep TIFFTAG_DATETIME
2
  • 1
    +1 for PostGIS. With PostgreSQL you can also take advantage of the ARRAY data type to store the array of dates of arbitrary length. Commented Dec 7, 2015 at 20:16
  • intersting idea @alpha-beta-soup ... not sure how qgis handles array types from postgres layers, but OP could always use this answer to 'explode' the array into one row per point/raster match (not quite sure what the proper term for that is...)
    – Steven Kay
    Commented Dec 7, 2015 at 21:54
3

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()
3

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.