0

I have multiple (381) high resolution GeoTIFF's (2x2 meter resolution, 1 band, 700 MB each) that I somehow want to upload to a PostGIS-server to be able to run ST_intersects or ST_3Dintersects queries to. I need the resolution preserved since every 2x2 meter square holds a height value in meters that new buildings can have before they are visible from a number of cultural heritage places. And to be extra clear the desired intersect operation I want to run is in the z-dimension, not only based on the coordinates. Ideally I would like to be able to create a polygon shapefile with a z-attribute for a planned building and run a st_intersect query to see from which of the 381 points the new building would be visible from.

When asking about this at Stackexchange before, the question got closed because it was to broad(which it indeed was). One user questioned the need to load the raster's to postgis when the footprint of the image was all I needed. Since then I've been trying to understand what a raster footprint is, but I'm not sure I've figured it out yet:

I've tried to Polygonize a raster with QGIS (Raster>Conversion>Polygonize) which worked and stores the z-values I want to use BUT the vector version of the 700 MB GeoTIFF takes up 4,33 GB of storage. So I'm guessing this approach is even worse then uploading the raster straight to PostGIS with raster2pgsql? So I'm guessing this is not what the user had in mind when suggesting footprints as a way forward?

I've also tried a plugin for QGIS called Image Footprint which from the plugin description sound like just the plugin I need:

This plugin create a catalog layer from directories of images, how a memory layer, where each feature has metadata and footprint of a image. See presentation in http://pt.slideshare.net/LuizMotta3/qgis-ibama-imagefootprint

The plugin gives me two options: 1. Bounding Box or, 2. Valid Pixels.

  1. The bounding box shapefile is only 176 KB big which is handy, has two columns "filename" and "meta_html", but no values from the raster is present in the shapefile which only contains one feature and some meta data about the area of valid pixels, number of holes, projection etc.

  2. The Valid pixels shapefile is also 176 KB and has almost the same meta data except for more info about pixel size and extent.

But none of the above results from this plugin gives me any actual info from the 2x2 meter cells in the original raster. Although the attributes contain a path to where the actual raster is stored on the hard drive. So this can not be the way forward either?

My question is how would you store these high resolution GeoTIFF's on a PostGIS-server to be able to run the above mentioned st_intersects query with the fastest possible performance?

Server runs Ubuntu 17.10 with Postgresql 9.6.6/PostGIS 2.4.1 on a Intel i5-6600k CPU with 4 cores and 4 logical processors, 16 GB om physical RAM and a 1 TB SSD harddrive.

Metadata from one of the rasters: enter image description here

Improve this question
4
  • 1
    You still have multiple questions here. A footprint is the quadrilateral polygon constructed from corner points e.g. { SWx, SWy, SEx, SEy, NEx, NEy, NWx, NWy, SWx, SWy}). Loading the rasters into the database solely for the purpose of determining overlap is like using an A380 to deliver mail in the bush.
    – Vince
    Commented Jan 8, 2018 at 17:57
  • Thanks for the first explanation I've received on what a footprint is, very helpful! I removed the duplicate question about raster2pgsql. If I understand your definition of a footprint correctly, loading a footprint won't do it for me since I need all the z-values in every 2x2 meter cell of the 700 MB raster!? Looks like I'm gonna need an A380 to deliver the mail?
    – johlund
    Commented Jan 8, 2018 at 19:40
  • No, you missed it again. The footprint can identify which image, at which point any tool can identify the pixels.
    – Vince
    Commented Jan 8, 2018 at 20:04
  • Sorry Vince! So the path-column that the plugin Image Footprint creates in the footprint attribute table can be used by PostGIS to get me the answer without having to load a singel raster file into the db? Even though it's not solely a xy, but mainly a z-overlap I wanna determine?
    – johlund
    Commented Jan 8, 2018 at 20:15

1 Answer 1

Reset to default
1

I'm not sure why you've been told to upload the footprints.

If your elevation data does not overlap, you can upload all the rasters into one big table, and it will probably be fast enough.

To load your data into Postgres, you need to use raster2pgsql. See this documentation.

The first time you do this, you'll want to use the default -c flag, and then subsequent data should be loaded with -a to append.

It's a pretty big job, and there's a lot of tuning you can do. For example, when you store raster data in PostGIS, the loader tiles it up for you. I ended up using a 20x20 cell grid for the data that I loaded. And this way you have a row for each 20x20 'tile'. You create a spatial index on the bounds of each tile, and then you take this into account when you run a query. I.e., you select out the values at a lat/lon from a tile that intersects those lat/lon to get a single discrete value. Or you select all the tiles that intersect an area and merge them together into a TIFF to grab a subset out.

Improve this answer
5
  • Thank you so much Alex for a kind and very useful answer! My GeoTIFF's does overlap since they are viewshed results over the city of Stockholm, just from 381 different points. Does this make tiling unnecessary (it doesn't speed up the queries?) and would you still recommend uploading to the same table?
    – johlund
    Commented Jan 9, 2018 at 6:56
  • Right now my biggest concern is that there (up to this date) is no working intersect query in PostGIS that can handle z intersections based on a DEM/GeoTIFF. Just a TIN: postgis.net/docs/ST_3DIntersects.html
    – johlund
    Commented Jan 9, 2018 at 7:03
  • Hi @johlund, so your query will be for YES/NO for whether a new point is visible from one of the 381 different points? You could definitely form a query for this with all the data in one table. In fact, that should be the fastest way to do it.
    – Alex Leith
    Commented Jan 9, 2018 at 22:21
  • Alex, thank you so much for positive feedback! Your definition above is probably the most simple/straight-forward way to explain my problem. But in my case the defintion would be something a litte bit more complex like: YES/NO for whether a new polygon with z-values spatialy intersects with one of the 381 different viewshed rasters? And there is still the problem with the raster format. St_3Dintersects handles TIN's and my rasters are GeoTIFFs (300 GB in total so converting will be time consuming)
    – johlund
    Commented Jan 10, 2018 at 9:11
  • Do you need a 3D instersection? Maybe you can prepare a query that finds the raster cells that intersect a polygon, and then compare the max or min heights in the cells to the polygon? I think you'd better start with a proof of concept, so just load one raster and see if you can come up with a query that works.
    – Alex Leith
    Commented Jan 10, 2018 at 22:54

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