4

I have a set of point geometries that I would like to cluster together using ST_ClusterDBSCAN. Unfortunately, it doesn't seem to work for WGS84 and requires a meter-based projection. Every tutorial that I can seem to find manually determines which projection is a good meter-based projection for their particular dataset. Is there a generalizable way to do this automatically?

Basically what I'd like to do is take the convex hull of my points, take the centroid of that polygon and then query to find what's a good meter-based projection for that particular region of the world. I'll then convert my points to this projection and run the clustering.

Matt
  • 402
  • 3
  • 10
  • What is the 'bounding box' of your data? That is, what area does your data occupy - that will help you determine a UTM zone, for example - which are in metres. Here are the zones: http://www.dmap.co.uk/utmworld.htm – DPSSpatial_BoycottingGISSE Mar 16 '18 at 14:58
  • OK, so basically I want to get the UTM zone for a given point. If i can figure out which zone the point belongs to like on this grid, then how would I get the appropriate SRID? – Matt Mar 16 '18 at 15:36
  • If your data spans more than 1 UTM zone, you might want to check into something like Albers Equal Area Conic, which I believe is in Metres... pinging each UTM zone might be tricky... – DPSSpatial_BoycottingGISSE Mar 16 '18 at 15:49
  • you can probably rework the answer from https://gis.stackexchange.com/questions/269518/auto-select-suitable-utm-zone-based-on-grid-intersection – Ian Turton Mar 16 '18 at 16:01
  • also see http://cite.opengeospatial.org/OGCTestData/wms/1.1.1/spec/wms1.1.1.html#auto_projections.42001 for a way of auto choosing a UTM projection – Ian Turton Mar 16 '18 at 16:07

2 Answers2

4

If you're willing to wander into the murky world of undocumented PostGIS functions, you could consider _ST_BestSRID. It accepts a geography object and spits back a recommended SRID. The code can be seen here and is well-commented (scroll ahead to line 843) so you can get an idea of the logic used.

Note that the SRIDs for UTM zones returned by this function are internal to PostGIS. They're not EPSG codes and won't appear in spatial_ref_sys but are defined in source here.

dbaston
  • 13,048
  • 3
  • 49
  • 81
  • Thanks @dbaston, I came across ST_BestSRID here and was hoping to avoid it, but if there is no other option, I guess that'll have to do. – Matt Mar 16 '18 at 17:33
  • Do you prefer to avoid it because of the lack of documentation, or are you looking for different logic to select the SRID? – dbaston Mar 16 '18 at 17:49
  • I just don't want to update my version of PostGIS and have it not be there or be called something else. I guess I just figured it's private and not-documented for good reason. – Matt Mar 16 '18 at 18:10
  • @dbaston the here link is broken would you mind updating it? I am trying to find out what these codes mean – ziggy Nov 13 '19 at 17:11
1

Update based on Ian's comment below...

Use US National Atlas Equal Area, assuming your data is in North America

http://spatialreference.org/ref/epsg/2163/

This coordinate system is already available in PostGIS.

And here are the specs that detail the coordinate system units in metres:

http://spatialreference.org/ref/epsg/2163/html/

PROJCS["US National Atlas Equal Area",
    GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",
        DATUM["Not_specified_based_on_Clarke_1866_Authalic_Sphere",
            SPHEROID["Clarke 1866 Authalic Sphere",6370997,0,
                AUTHORITY["EPSG","7052"]],
            AUTHORITY["EPSG","6052"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4052"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",45],
    PARAMETER["longitude_of_center",-100],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","2163"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]
DPSSpatial_BoycottingGISSE
  • 18,790
  • 4
  • 66
  • 110