6

Question Summary

I am trying to find all the MultiPolygons within lat/lng bounds, and I am unable to get a valid query. Ultimately I believe I have not transformed the MultiPolygons into something i can query correctly, but I might be reaching...

Question Details

NOTE: I have provided all of the Background Details of my data at the end of the post, and demonstrated (hopefully) that it is correct and has the appropriate SRIDs .

I have two GEOMETRY columns:

geom        GEOMETRY(MultiPolygon,2926) => original from the import
geom_google GEOMETRY(MultiPolygon,4326) => created, so i could query

Assuming i have the correct data and SRIDs, I am trying to get all the rows that are within a certain bounds

The bounds are coming from a google map getBounds() call, and are:

{
    north: 48.04148960598052,   // lat
    south: 47.18983904097505,   // lat
    east: -121.80370233688086,  // lng
    west: -122.65239618453711   // lng
}

To simplify:

{
    north: 49,   // lat
    south: 47,   // lat
    east: -121,  // lng
    west: -123   // lng
}

Points are set as ST_SetSRID(ST_Point(lng, lat),SRID)

ST_MakeBox2D asks for: pointLowLeft and pointUpRight, so:

pointLowLeft => west(lng) and south(lat)
    ST_SetSRID(ST_Point(-123, 47),4326)

pointUpRight => east(lng) and north(lat): 
    ST_SetSRID(ST_Point(-121, 49),4326)

My Box would then be:

ST_MakeBox2D(
    ST_SetSRID(ST_Point(-123, 47), 4326),
    ST_SetSRID(ST_Point(-121, 49), 4326)
)

which when selected shows:

BOX(-123 47,-121 49)

Attempts at an Answer

select * from "kc"."plss" where ST_Intersects(
    ST_MakeBox2D(
        ST_SetSRID(ST_Point(-123, 47), 4326),
        ST_SetSRID(ST_Point(-121, 49), 4326)
    ),
    geom_google);

which returns:

[XX000] ERROR: Operation on mixed SRID geometries

I have tried:

select * from "kc"."plss" where ST_Intersects(
  ST_Transform(
      ST_MakeBox2D(
        ST_SetSRID(ST_Point(-123, 47), 4326),
        ST_SetSRID(ST_Point(-121, 49), 4326)
      ),
      4326
  ),
  ST_SetSRID(geom_google, 4326));

which returns:

[XX000] ERROR: Input geometry has unknown (0) SRID

I have also tried:

select ST_AsGeoJSON((ST_DUMP(ST_Transform(geom_google,4326))).geom::geometry(Polygon,4326)) from "kc"."plss" where ST_Intersects(ST_MakeBox2D(
    ST_SetSRID(ST_Point(-123, 47), 4326),
    ST_SetSRID(ST_Point(-121, 49), 4326)
  ),geom);

Which is valid, but returns 0 rows

Actual Question

Can anyone point me to where I am going wrong?


FINAL, CORRECT QUERY:

Here is the final, correct query (helped by @tilt), just so it is clear:

SELECT id from "kc"."plss"
WHERE ST_Intersects(
ST_SetSRID(
    ST_MakeBox2D(
        ST_Point(-122.23124645385474, 47.620059604358204),
        ST_Point(-122.20472477111548, 47.64666550757364)
    ),
    4326
),
geom_google);

which results in 20 ids:

1090,1096,1093,1095,1178,1179,1180,1177,1231,1234,1235,1236,1313,1314,1315,1316,1367,1369,1364,1366

or for those of you looking for the polygons, here is my real query:

SELECT ST_AsGeoJSON((ST_DUMP(geom_google)).geom::geometry(Polygon,4326)) from "kc"."plss"
WHERE ST_Intersects(
    ST_SetSRID(
        ST_MakeBox2D(
            ST_Point(-122.23124645385474, 47.620059604358204),
            ST_Point(-122.20472477111548, 47.64666550757364)
        ),
        4326
    ),
    geom_google);

Background Details

The dataset I imported is:

http://www5.kingcounty.gov/sdc/Metadata.aspx?Layer=plss

This data is in http://spatialreference.org/ref/epsg/2926/

I imported it into my db, and created a new column:

ALTER TABLE kc.plss ADD geom_google GEOMETRY(MultiPolygon,4326) NULL;

Then ran:

UPDATE "kc"."plss"
SET  geom_google = ST_Transform(geom,4326);

Now if I query:

select ST_SRID(geom), ST_SRID(geom_google) from "kc"."plss" LIMIT 1;

I get:

2926    4326

And then if I do:

select ST_AsGeoJSON(geom), ST_AsGeoJSON(geom_google) from "kc"."plss" LIMIT 1;

I get (truncated, of course):

{"type":"MultiPolygon","coordinates":[[[[1246963.803,247960.583000004],[1246920.789,246646.928000003],...]]]]}
{"type":"MultiPolygon","coordinates":[[[[-122.430314129407,47.6685817588591],[-122.430378032134,47.6649787218942],...]]]]}

So, things are looking good.

I also added some indexes:

CREATE INDEX ix_spatial_geom ON "kc"."plss" using gist(geom);
CREATE INDEX ix_spatial_geom_google ON "kc"."plss" using gist(geom_google);
  • Compliments on the nicely elaborated question. – tilt Jun 1 '17 at 7:44
  • I accept that compliment! ;-) This isn't my first rodeo, and I hoped to get it answered. Nothing beats a well phrased question to get an answer! Thank you for helping. – Joshua Ball Jun 1 '17 at 15:17
  • for the inner point constructors, you can remove ST_SetSRID calls, and just use ST_Point like the example in the manual – Mike T Jun 1 '17 at 22:56
  • your question should be used as a template on how to properly ask a sound question, with edits, follow ups and good structure – ziggy Jun 2 '17 at 17:47
4

The problem is that ST_MakeBox2D erases the srid. So you have to do st_setsrid after you did makebox2d (also see the example in the postgis docs: https://postgis.net/docs/ST_MakeBox2D.html

  • Indeed you are correct. In fact, I had at one time wrapped ST_MakeBox2D in ST_SetSRID in one the combinations of the queries I had written, but at the time I believe I was failing for another reason (my box was mixing up the points lowLeft and upRight maybe?), and so I had failed for other reasons. – Joshua Ball Jun 1 '17 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.