postgis shape file import problems

Question

Hi I'm trying to import a shape file from

http://www.nyc.gov/html/dcp/html/bytes/bytesarchive.shtml

into a postgis database. the above files creates MULTIPOLYGONS when i import using shp2pgsql.

then i'm trying to simply determine if lat/long points are contained in my multipolygons

however my select's are not working, and when i print out the poitns of my the_geom column it seems to be very broken.

select st_astext(geom) from (select (st_dumppoints(the_geom)).* from nybb where borocode =1) foo;

gives the result...

            st_astext
------------------------------------------
 POINT(1007193.83859999 257820.786899999)
 POINT(1007209.40620001 257829.435100004)
 POINT(1007244.8654 257833.326199993)
 POINT(1007283.3496 257839.812399998)
 POINT(1007299.3502 257851.488900006)
 POINT(1007320.1081 257869.218500003)
 POINT(1007356.64669999 257891.055800006)
 POINT(1007385.6197 257901.432999998)
 POINT(1007421.94509999 257894.084000006)
 POINT(1007516.85959999 257890.406100005)
 POINT(1007582.59110001 257884.7861)
 POINT(1007639.02150001 257877.217199996)
 POINT(1007701.29170001 257872.893099993)
...

for points in nyc, this is very off.. what am i doing wrong?

This question will probably fit better on http://gis.stackexchange.com You might consider asking a moderator for it to be moved there. — steenhulthin, Aug 21 '11 at 12:40

score 4 · Accepted Answer · answered Aug 21 '11 at 13:01

The points are not of. The spatial data that is referred to is NOT in lat/long. This is why numbers are different from what you expect. If you need it to be in long/lat it must be reprojected. See more here: http://postgis.refractions.net/news/20020108/

The projection of the data seems to be in the NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet coordinate system (according to the metadata - see code.).

<spref>
    <horizsys>
        <planar>
            <planci>
                <plance Sync="TRUE">coordinate pair</plance>
                <coordrep>
                    <absres Sync="TRUE">0.000000</absres>
                    <ordres Sync="TRUE">0.000000</ordres>
                </coordrep>
                <plandu Sync="TRUE">survey feet</plandu>
            </planci>
            <mapproj><mapprojn Sync="TRUE">Lambert Conformal Conic</mapprojn><lambertc><stdparll Sync="TRUE">40.666667</stdparll><stdparll Sync="TRUE">41.033333</stdparll><longcm Sync="TRUE">-74.000000</longcm><latprjo Sync="TRUE">40.166667</latprjo><feast Sync="TRUE">984250.000000</feast><fnorth Sync="TRUE">0.000000</fnorth></lambertc></mapproj></planar>
        <geodetic>
            <horizdn Sync="TRUE">North American Datum of 1983</horizdn>
            <ellips Sync="TRUE">Geodetic Reference System 80</ellips>
            <semiaxis Sync="TRUE">6378137.000000</semiaxis>
            <denflat Sync="TRUE">298.257222</denflat>
        </geodetic>
        <cordsysn>
            <geogcsn Sync="TRUE">GCS_North_American_1983</geogcsn>
            <projcsn Sync="TRUE">NAD_1983_StatePlane_New_York_Long_Island_FIPS_3104_Feet</projcsn>
        </cordsysn>
    </horizsys>
</spref>

If you work much with spatial data I suggest that you read more about map projection.

Use the -s flag in shp2pgsql to set the SRID to [2263](http://prj2epsg.org/epsg/2263). Then you can use [ST_Transform](http://postgis.refractions.net/docs/ST_Transform.html) to convert to 4326. — Paul Ramsey, Aug 22 '11 at 18:04
Another possibility - pass ```-s 2263:4326``` to shp2pgsql. Does the conversion for you. — Alastair, Jun 02 '13 at 04:24

Grzegorz Szpetkowski · Answer 2 · 2011-08-21T12:32:20.670

I think this is not issue with PostGIS. I checked input esri Shape file nybb.shp with AvisMap Free Viewer and as you see points are weird itself:

enter image description here

However there is something interesting in nybb.shp.xml metadata file:

<spdom>
    <bounding>
        <westbc Sync="TRUE">-74.257465</westbc>
        <eastbc Sync="TRUE">-73.699450</eastbc>
        <northbc Sync="TRUE">40.915808</northbc>
        <southbc Sync="TRUE">40.495805</southbc>
    </bounding>
    <lboundng>
        <leftbc Sync="TRUE">913090.770096</leftbc>
        <rightbc Sync="TRUE">1067317.219904</rightbc>
        <bottombc Sync="TRUE">120053.526313</bottombc>
        <topbc Sync="TRUE">272932.050103</topbc>
    </lboundng>
</spdom>

I am not familiar with those toolkit (ESRI ArcCatalog), but most probably you need to rescale your points after import using that metadata.

postgis shape file import problems

2 Answers2