0

I'm trying to ST_Clip() raster in `SRID:32617' with a polygon.
What I did

  • ST_Transform(polygon, 32617): It fails, couldn't project point() exceeded limit (-14)
  • ST_Transform(rast, 4326, 'NearestNeighbor'): It runs really slow (more than 700,000,000 ms, then I terminated before it finished) about 1300 rasters (tiles)

My polygon got lon/lat from Google map:

ST_GeomFromText('POLYGON((105.78873 21.38659, 106.49185 21.38659, 105.78873 20.71758, 
106.49185 20.71758, 105.78873 21.38659))', 4326)

I'm using postGIS 2.1

UPDATE:
My raster is retrieved from Landsat8. I downloaded the raster from earthexplorer.usgs.gov
The raster is imported into postGIS by raster2pgsql with tile-size 500x500, SRID is kept as original (32617,32617, 32632, 32648, 32649, 32658, 32659, 32660))

UPDATE:
Here's Map projection in metadata text of one scene.

GROUP = PROJECTION_PARAMETERS
MAP_PROJECTION = "UTM"
DATUM = "WGS84"
ELLIPSOID = "WGS84"
UTM_ZONE = 49
GRID_CELL_SIZE_PANCHROMATIC = 15.00
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
ORIENTATION = "NORTH_UP"
RESAMPLING_OPTION = "CUBIC_CONVOLUTION"
END_GROUP = PROJECTION_PARAMETERS

UPDATE The polygon above is invalid. Here's correction: ST_GeomFromText('POLYGON((105.78873 21.38659, 106.49185 21.38659, 106.49185 20.71758, 105.78873 20.71758, 105.78873 21.38659))', 4326)

10
  • I deleted my original answer. Please update your question with more details. What is the source of your raster? Where did you download it? What is the spatial reference system information for it? How did you import it into postgis?
    – BradHards
    Commented Mar 11, 2014 at 6:36
  • Hi @BradHards, I updated my question.
    – Jackie
    Commented Mar 11, 2014 at 6:50
  • 1
    Why do you think your scene is in 32617? The metadata you show says UTM Zone 49.
    – BradHards
    Commented Mar 11, 2014 at 7:10
  • After importing it into PostGIS, I run ST_SRID(rast) and it return 7 SRIDs (32617, 32632, 32648, 32649, 32658, 32659, 32660). Because, I have about 60 scenes, some of them have SRID:32617
    – Jackie
    Commented Mar 11, 2014 at 7:31
  • 1
    How can an image in 32617 match a polygon at 106E?
    – BradHards
    Commented Mar 12, 2014 at 5:23

1 Answer 1

2

That spatial reference identifier (32617) is in UTM Zone 17N (which is in the US). It makes no sense at all to try to project something at 106E to Zone 17N (78W to 84W) coordinates.

You probably want to use something like UTM Zone 48N: http://spatialreference.org/ref/epsg/32648/

SELECT AsEwkt(ST_Transform(ST_GeomFromText('POLYGON((105.78873 21.38659, 106.49185 21.38659, 105.78873 20.71758, 106.49185 20.71758, 105.78873 21.38659))', 4326), 32648));
SRID=32648;POLYGON((581760.7532212715 2365138.811852233,654656.3711929469 2365667.925132345,582126.9328327369 2291091.94133294,655349.2256348798 2291607.543684119,581760.7532212715 2365138.811852233))

[I used spatialite, but postgis will likely look similar]

Now you might have rasters in another UTM zone, but you can probably assume that if the srid doesn't match, then the raster is in another UTM zone and it won't match any polygon in this UTM zone.

3
  • Hi @BradHards, I'm trying to clip the raster with a polygon, so I need to make them share the same SRID. However, ST_Transform(rast, 32648, 'NearestNeighbor') takes extremely long time. Any suggestion? Thanks for your replying
    – Jackie
    Commented Mar 11, 2014 at 4:18
  • This is a different question. I think I've answered the question you asked.
    – BradHards
    Commented Mar 11, 2014 at 6:17
  • Thanks for your time, but I still couldn't ST_Clip() raster SRID:32617 with a polygon - this is my question.
    – Jackie
    Commented Mar 11, 2014 at 6:33

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.