3

I have done:

  • Import a SHP file into PostGIS by PostGIS SHP loader GUI.

  • Import a raster into PostGIS:

    \raster2pgsql -I -s 4326 landsat.TIF -t 500x500 public.landsat > query.sql

After importing all into postGIS, I have tried exporting raster to PNG by ST_AsPNG. It run well. Besides, I run the example of ST_Clip() in PostGIS docs. It worked.

Question: I am trying to clip a raster with a polygon of the SHP file. It returns nothing.
I tried ST_Intersects(rast, geom polygon) like this. It failed.

Can anyone help me on clipping a raster with a pylygon of SHP file?

Viewing in QGIS, the raster overlaps the polygon. I'm using PostGIS 2.0.

EDIT: Here's what I do with ST_Intersects()

SELECT 
ST_Intersects(polygon.geom, raster.rast)
FROM 
  public.vnm_adm2 as polygon,
  public.landsat as raster;

And it returns false in all rows.

UPDATE: There're some things to be noted here:

  • I used ST_AsPNG(), which remove spatial reference to output file.
  • Switch -s [4326] doesn't work yet (PostGIS 2.0), so the rasters ain't be converted to polygon's EPSG, then it returns false in all rows.
Community
  • 1
Jackie
  • 303
  • 1
  • 3
  • 10
  • When you say "It failed", can you add some more detail on the exact calls that you used (from the postgis documentation example, and from your shapefile example) , and the exact result for the shapefile example? That is, we know it didn't work, but we don't know the error message. Just click edit below the question to add this critical detail. – BradHards Jan 24 '14 at 03:17
  • I used the ST_Clip() example from PostGIS 2.0 document. I have updated the question. – Jackie Jan 24 '14 at 03:32

1 Answers1

4

Your ST_Intersects didn't failed, because - as you said - it returns false in all rows!

It fails when you get an error, but this isn't your case!

So, why did you get false, even your layers overlap in QGIS? Well, this is a tricky one...

What you see in QGIS, may be a polygon layer loaded from a PostGIS table, but I'm pretty sure that your raster is loaded from a local folder!

Please, try this exercise in QGIS:

  • right click on the polygon layer, then choose Properties;
  • in the General tab, be sure to have your PostGIS database table as a layer source.
  • be sure your layer Coordinate Reference System is EPSG:4326.

Repeat these steps for the raster layer.

If you didn't see anything wrong then you may skipped something when you loaded the raster.

You'll have to figure out yourself where you went wrong, and, to help you easier identify the error, I've prepared for you a fully functional project. Inside the project zip, in a readme file, you'll find all the necessary steps to obtain true, after a polygon with raster ST_Intersects operation.

There is a georeferenced tiff raster, which can be loaded in your PostGIS:

raster& polygon.png

The raster is splitted (through the raster2pgsql command) in 3 x 4 tiles and stored in a raster table using 12 rows. If you'll follow the instructions you'll see that everything works as expected:

first_intersects.png

This is how you find if the polygon intersects a specific raster tile:

enter image description here

Sorin Călinică
  • 5,118
  • 1
  • 17
  • 23
  • Hi @Sorin, sorry about late replying. New Year holiday has just end here. Thanks for your example project, it helps me to understand PostGIS a lot. I followed your instruction and got the same result as you said.
    Besides, I tried ST_Clip() and I got what I expected.
    I will apply in my project and talk to you later.
    Thank again :)     – Jackie Feb 10 '14 at 04:08
  • You're right @Sorin. My problem is that I tried to transform raster original georeference to other when using raster2pgsql. For detail, my raster is in 32649, but I used: raster2pgsql -s 4326. It leads to the result the ST_Intersects() returns false in all rows. Thank you so much! – Jackie Feb 10 '14 at 08:24
  • 1
    Should I update my question to show where I went wrong and how to solve it! – Jackie Feb 10 '14 at 09:10
  • 1
    I suggest you to post your solution as an answer, instead of editing the question. It would be more useful for other user with similar problems. – Sorin Călinică Feb 10 '14 at 09:15