[postgis-users] Question about using ST_Transform()

[Pierre Racine]

I generally display things in OpenJump to make sure they align correctly before doing anything instead of assuming they do. Do the projected geometries align properly with the rasters?

You can display the raster extents easily in OpenJump by casting the raster to geometry (rast::geometry)...

You can also do this easily in QGIS or ArcGIS 10.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
> bounces at postgis.refractions.net] On Behalf Of David Quinn
> Sent: Thursday, February 09, 2012 10:45 AM
> To: PostGIS Users Discussion
> Subject: [postgis-users] Question about using ST_Transform()
> 
> Hello All,
> 
> I'm trying to understand how to use ST_Transform() on a query. When I have
> everything in the same projection (in this case Albers Equal Area Conic) my query
> returns the correct result:
> 
> SELECT
>     SUM((stats).sum)
> FROM (
>     select
>         ST_SummaryStats(ST_Clip(r.rast, 1, p.geom, NULL, TRUE)) as stats
>     FROM
>         ny_albers as r,
>         cities_albers as p
>     WHERE
>         p.name = 'new_york'
>     AND
>         ST_Intersects(r.rast, p.geom)
>     )  as foo;
> 
> If I use a table with polygons in a different projection (in this case Web
> Mercator), I do not get any values returned:
> 
> SELECT
>     SUM((stats).sum)
> FROM (
>     select
>         ST_SummaryStats(ST_Clip(r.rast, 1, ST_Transform(p.geom,102003), NULL,
> TRUE)) as stats
>     FROM
>         ny_albers as r,
>         cities_webMercator as p
>     WHERE
>         p.name = 'new_york'
>     AND
>         ST_Intersects(r.rast, ST_Transform(p.geom,102003))
>     )  as foo;
> 
> I added the projection SRID:102003 using the following lines:
> 
> INSERT into
>     spatial_ref_sys (
>     srid,
>     auth_name,
>     auth_srid,
>     srtext,
>     proj4text
>     )
> values (
>      102003, -- I changed this line as it was originally 9102003 which did not satisfy
> the constraints
>      'esri',
>      102003,
> 'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_A
> merican_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980
> ",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453
> 292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_
> Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center
> ",-
> 96],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2
> ",45.5],PARAMETER["latitude_of_center",37.5],UNIT["Meter",1],AUTHORITY["E
> PSG","102003"]]',
>      '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0
> +ellps=GRS80 +datum=NAD83 +units=m +no_defs ');
> 
> I had some problems when I tried to add this projection from
> http://spatialreference.org/ref/esri/102003/postgis/
> The names and values of srtext, and proj4text were reversed, and the SRID
> seemed incorrect so I changed it. I don't know if this is related to my question
> above, but I thought it worth mentioning. I did import the raster and polygon
> using this SRID on the query that worked.
> 
> Thank you,
> David
>