0

tldr; I cannot get ST_Within to find business points that fall within two sets of polygons. I keep getting zero results. All three tables have geom col with SRID 4326.

I have a data table where I have two columns prem_city_ps (LA city boundaries) and prem_council (LA council districts ). I populate those columns using the two ST_Within statements below which I run in consecutive order.

My problem is I can only get 'UPDATE 0' results when running those ST_Within statements.

I know the lat,lon points in the data are accurate. I create a 'geom' column with:

UPDATE poursafe.ca_licenses SET geom = ST_SetSRID(ST_MakePoint(lon,lat),4326);

I've made sure all three datasets are SRID 4326, and confirm with:

SELECT Find_SRID('poursafe', 'ca_licenses', 'geom');

Here's sample data: https://github.com/PourSafe/data_samples including:

  1. Sample ca_licenses data
  2. Sample ca_la_la_council data Source: https://egis3.lacounty.gov/dataportal/2012/08/07/la-city-council-districts-2012/
  3. Sample ca_la_city_boundaries data Source: https://data.lacounty.gov/GIS-Data/City-Boundaries/wyja-ysz2

Here are my create table statements

Create Data Table

-- Name: ca_licenses; Type: TABLE; Schema: poursafe; Owner: postgres

CREATE TABLE poursafe.ca_licenses (
    id bigint DEFAULT nextval('poursafe.licenses_id_seq'::regclass) NOT NULL,
    license character varying NOT NULL,
    lat double precision,
    lon double precision,
    geom public.geometry(Point,4326),
    prem_city_ps character varying,
    prem_council character varying
);

ALTER TABLE poursafe.ca_licenses OWNER TO postgres;

I created an index of the data's geom col with:
CREATE INDEX ca_licenses_gix ON poursafe.ca_licenses USING GIST (geom);

Other tables were created in the PostGis Shapefile Import/Export Manager process indicated below.

QUERIES

Note: I manually created indexes for these two public tables because I'm not sure how to test the ones created by PostGis Shapefile Import/Export Manager (and dropped those created by that tool).

The first adds a '49' in prem_city_ps to any record who's geom point falls within these many polygons). This is the dataset's ID number for LA:

1. Upload shape file from: https://data.lacounty.gov/GIS-Data/City-Boundaries/wyja-ysz2
    SRID 4326, options: create spatial after load, load via 'copy'
2. Set SRID to 4326:
    SELECT UpdateGeometrySRID('public','ca_la_city_boundaries','geom',4326);
3. SELECT Find_SRID('public', 'ca_la_city_boundaries', 'geom');
4. CREATE INDEX ca_la_city_boundaries_gix ON public.ca_la_city_boundaries USING GIST (geom);

Update prem_city_ps col:

UPDATE poursafe.ca_licenses dst
SET prem_city_ps = 'LOS ANGELES'
WHERE EXISTS (
  SELECT src.geom, src.city_no
  FROM public.ca_la_city_boundaries src
  WHERE ST_Within(dst.geom, src.geom )
  AND src.city_no = '49'
  );

The second enters the district number in the prem_council col for records who's geom point falls within these 15 city council polygons:

1. Upload shape file from: https://egis3.lacounty.gov/dataportal/2012/08/07/la-city-council-districts-2012/
    SRID 2229, options: create spatial after load, load via 'copy'

[[PostGis Shapefile Import/Export Manager did not set SRID in upload]] 2. Set SRID to 2229: SELECT UpdateGeometrySRID('public','ca_la_la_council','geom',4326); 3. change from SRID 2229 to 4326: ALTER TABLE public.ca_la_la_council ALTER COLUMN geom TYPE geometry(MultiPolygon,4326) USING ST_Transform(ST_SetSRID(geom,2229),4326); 4. Check SRID with: SELECT Find_SRID('public', 'ca_la_la_council', 'geom'); 5. CREATE INDEX ca_la_la_council_gix ON public.ca_la_la_council USING GIST (geom);

update (prem_council) col:

UPDATE poursafe.ca_licenses dst
    SET prem_council = district
FROM public.ca_la_la_council src
WHERE ST_Within(dst.geom,src.geom) AND dst.geom&&src.geom;

EXPLAIN ANALYZE RESULTS

EXPLAIN ANALYZE
UPDATE poursafe.ca_licenses dst
    SET prem_council = district
FROM public.ca_la_la_council src
WHERE ST_Within(dst.geom,src.geom) AND dst.geom&&src.geom;

"Update on ca_licenses dst  (cost=0.13..68492.30 rows=1 width=390) (actual time=1351.828..1351.828 rows=0 loops=1)"
"  ->  Nested Loop  (cost=0.13..68492.30 rows=1 width=390) (actual time=0.316..1250.973 rows=9015 loops=1)"
"        ->  Seq Scan on ca_licenses dst  (cost=0.00..19858.80 rows=117880 width=316) (actual time=0.012..124.320 rows=117880 loops=1)"
"        ->  Index Scan using ca_la_la_council_gix on ca_la_la_council src  (cost=0.13..0.40 rows=1 width=106) (actual time=0.008..0.009 rows=0 loops=117880)"
"              Index Cond: ((geom ~ dst.geom) AND (dst.geom && geom))"
"              Filter: _st_contains(geom, dst.geom)"
"              Rows Removed by Filter: 0"
"Planning time: 0.980 ms"
"Execution time: 1351.911 ms"


EXPLAIN ANALYZE
UPDATE poursafe.ca_licenses dst
SET prem_city_ps = 'LOS ANGELES'
WHERE EXISTS (
  SELECT src.geom, src.city_no
  FROM public.ca_la_city_boundaries src
  WHERE ST_Within(dst.geom, src.geom )
  AND src.city_no = '49'
  );

"Update on ca_licenses dst  (cost=0.14..69626.74 rows=39 width=354) (actual time=5711.997..5711.997 rows=0 loops=1)"
"  ->  Nested Loop Semi Join  (cost=0.14..69626.74 rows=39 width=354) (actual time=2.391..5587.446 rows=9017 loops=1)"
"        ->  Seq Scan on ca_licenses dst  (cost=0.00..19858.80 rows=117880 width=316) (actual time=0.033..126.802 rows=117880 loops=1)"
"        ->  Index Scan using ca_la_city_boundaries_gix on ca_la_city_boundaries src  (cost=0.14..0.41 rows=1 width=38) (actual time=0.046..0.046 rows=0 loops=117880)"
"              Index Cond: (geom ~ dst.geom)"
"              Filter: ((city_no = '49'::numeric) AND _st_contains(geom, dst.geom))"
"              Rows Removed by Filter: 0"
"Planning time: 0.367 ms"
"Execution time: 5712.572 ms"

I (may clearly have already revealed) I'm relatively new to this. I have had this working before (see https://poursafe.com/los-angeles-city/ but, had to somewhat start over with my data.

  • 1
    the ca_la_la_council data, as per the linked source (not your GitHub repo), is projected in (or better, very close to) EPSG:2229...did you transform into EPSG:4326 in a previous step? if not; simply setting the CRS will render your data useless! btw., your repo is missing the sample ca_licenses data...,) – ThingumaBob Sep 9 '18 at 10:35
  • 1
    more recommendations: you should have (spatial) indexes in place on all tables...only then will the magic happen. the first UPDATE can easily be expressed like the second (the EXISTS construct might not even work as expected here), while in the second, the bbox comparison should probably come first (it should actually be unnecessary, since ST_Within is ought to include that already...you can find out if it does use the index by execute that query with EXPLAIN ANALYZE <query> and examine the plan...post it here also, maybe) – ThingumaBob Sep 9 '18 at 10:46
  • I uploaded the sample ca_licenses dataset in the github repo. I also added the queries section in my question which shows how I updated the sequence of uploading the shape file and setting/transforming the respective SRID's in my question. After doing that the queries worked like a charm. I posted the EXPLAIN ANALYZE results at the bottom of the question. I'm not sure how to promote your answer as the correct one as it's in comments. Thanks to everyone for the help. – enjoypb Sep 9 '18 at 23:42
  • well, I'd say, self answer your question; transform your updates into an answer (edit out of your question and add as an answer). that would make a little clearer what was going on. and you can even accept your answer, after 24h. – ThingumaBob Sep 13 '18 at 8:24
0

The order of parameters is not correct on the parameter list of ST_Within.

Extract from the PostGIS docs:

boolean ST_Within(geometry A, geometry B); 
Returns TRUE if geometry A is completely inside geometry B.

Swap the parameters:

ST_Within(src.geom, dst.geom)

or use the ST_Contains function

boolean ST_Contains(geometry geomA, geometry geomB);
Geometry A contains Geometry B if and only if no points of B lie in the exterior of A

Change the function name:

ST_Contains(dst.geom, src.geom)

  • I think dst are the points and src the polygons, in both UPDATEs; A inside B, points inside polygon...should be correct – ThingumaBob Sep 9 '18 at 10:37

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.