[postgis-users] error in geometry operations

[P Kishor]

On 8/20/06, Michael Fuhr <mike at fuhr.org> wrote:
> On Sun, Aug 20, 2006 at 09:06:08PM -0500, P Kishor wrote:
> > > The points table churned away for an hour and a half and
> > > successfully seems to have Transformed... but who knows,
> > > maybe it dint.
> >
> > # SELECT Extent(the_geom) FROM polys
> > "BOX(-6334962 -1864800.25,3377040.25 1407308.125)"
> >
> > # SELECT Extent(the_geom) FROM points
> > "BOX(-9914299 4535597.5,-4712810 13894887)"
> >
> > It sure dint. Both should be about the same, and both should cover the
> > US, Alaska, HI, and Puerto Rico. The polys seem correct, but coord
> > pairs for the points extent is whacky.
>
> Are you sure you committed the transaction that did the Transform?

As you know, I am using PgAdmin, and it is set for Autocommit as far
as I can tell (I don't have to commit anything, it gets committed).

> If you look at the points themselves (as opposed to using Extent),
> do they look transformed?

Now, that is the 5 million $ question. I have not looked at them
because I really have no good way of looking at them. I have uDig,
which mostly works on other datasets, but when I try to look at these
points, it croaks on me saying "Coordinate reference system unknown".
So, truthfully, I have no idea if they are transformed or not, let
alone transformed correctly. Here are the steps I followed to come
this far --

-- imported the lon/lat values into a table
COPY public.my_points (point_id,foo,longitude,latitude)
FROM 'C:/Documents and Settings/pkishor/My Documents/Data/my_points.txt'
WITH CSV HEADER

-- Point-ified the table
SELECT AddGeometryColumn('public', 'my_points', 'the_geom', 4326, 'POINT', 2);

-- made points from lon/lat (remembering that lon is 'x', lat is 'y')
UPDATE my_points
SET the_geom = SETSRID(MAKEPOINT(longitude, latitude), 4326);

-- indexed the geometry, and the ids for good measure
CREATE INDEX idx_myp_the_geom ON my_points USING GiST (the_geom);
CREATE INDEX idx_myp_point_id ON my_points USING btree (point_id);

-- wasted a few hours discovering the Proj.4 representation of srtext
-- for ESRI's projection definition. Added it to the spatial_ref table
INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text)
VALUES (
  102005,
  'ESRI',
  102005,
  'PROJCS[
    ''USA_Contiguous_Lambert_Conformal_Conic'',
    GEOGCS[
      ''GCS_North_American_1983'',
      DATUM[
        ''D_North_American_1983'',
        SPHEROID[''GRS_1980'',6378137.0,298.257222101]
      ],
      PRIMEM[''Greenwich'',0.0],
      UNIT[''Degree'',0.0174532925199433]
    ],
    PROJECTION[''Lambert_Conformal_Conic''],
    PARAMETER[''False_Easting'',0.0],
    PARAMETER[''False_Northing'',0.0],
    PARAMETER[''Central_Meridian'',-96.0],
    PARAMETER[''Standard_Parallel_1'',33.0],
    PARAMETER[''Standard_Parallel_2'',45.0],
    PARAMETER[''Latitude_Of_Origin'',39.0],
    UNIT[''Meter'',1.0]
  ];-19071146.8131673 -17820932.1267064 61.0351561931566;0 100000;0 100000;'
  '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0
+ellps=GRS80 +datum=NAD83 +units=m'
);

-- Changed the point srid to the new srs from the original 4326 lon/lat crap
UPDATE geometry_columns
SET srid = 102005 WHERE f_table_name = 'my_points';

-- mogrified the points permanently to the new projection
UPDATE my_points SET the_geom=Transform(the_geom, 102005);

-- discovered that it is all out of whack. Ran for help.
SELECT Extent(the_geom) FROM my_points
# "BOX(-9914299 4535597.5,-4712810 13894887)"


-- 
Puneet Kishor http://punkish.eidesis.org/
Nelson Inst. for Env. Studies, UW-Madison http://www.ies.wisc.edu/
Open Source Geospatial Foundation https://edu.osgeo.org/