2

I have installed PostGIS on my machine and I am using shapefiles provided by the US Census Bureau for finding zips and counties of coordinates. I am just using the PostGIS provided method for loading data.

I am testing my query using this data: https://gist.github.com/erichurst/7882666

This is my query:

select
  main.latitude,
  main.longitude,
  c.cntyidfp,
  s.statefp,
  z.zcta5ce
from (
  select
    latitude::numeric(15,8) as latitude,
    longitude::numeric(15,8) as longitude,
    st_setsrid(st_makepoint(longitude::numeric(15,8), latitude::numeric(15,8)), 4269) as mypoint
  from temp.sample
) as main
left outer join tiger.county as c
  on st_intersects(c.the_geom, main.mypoint)
left outer join tiger.state as s
  on s.statefp = c.statefp
left outer join tiger.zcta5 as z
  on st_intersects(z.the_geom, main.mypoint)

Now, all of the test data and the production data I will be using are coordinates (single points), but I keep getting this error:

ERROR: Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported. HINT: Change argument 1: 'GEOMETRYCOLLECTION(LINESTRING(-71.255183 42.736524,-71.255147 42.736554),POLY...'

It may be triggered by other points, this was only an example.

Why would a linestring be created? I am only testing points. I would like to be able to run the query on all of the data at once, not try each record individually. Is there something I can do to extract a point? Am I missing something? I am new to PostGIS but I have done a bit of research on this issue and others and I cannot seem to get around this.

  • Check the datatype of your zipcodes, county, state. One of them likely contains a line. It can be caused by small sliver between various administrative units. The example you provide is on the border between 2 states. – JGH May 18 '18 at 13:37
  • I see, I only checked data types of the sample data and they were (no surpise) all points. I will do that. If it is the case that some are indeed collections, is st_intersects the proper function? – prisonmike May 18 '18 at 13:40
  • No. The doc says: Do not call with a GEOMETRYCOLLECTION as an argument for geometry version. – JGH May 18 '18 at 13:42
  • 3
    st_intersects(st_collectionextract(z.the_geom,3),main.mypoint) – ziggy May 18 '18 at 13:44
  • Thanks for answering so quickly. I will test these out and update my answer when I can. I think your suggestions will save me a lot of time. Cheers! – prisonmike May 18 '18 at 13:55
1

Answering my own question here:

The issue that I am experiencing is related to how the ZCTAs are stored in the PostGIS database when created using the TIGER shapefiles.

From what I saw when looking through the data, the ZCTAs are stored as MULTIPOLYGONs, POLYGONs, or GEOMETRYCOLLECTIONs. When trying to use st_intersects/within/contains, the same error was popping up because of the GEOMETRYCOLLECTION objects.

One very slow solution is to check the type of each ZCTA geometry when running the query and extract when necessary (the solution given by ziggy). Another is to run the query on ZCTAs that are not GEOMETRYCOLLECTIONs and run another using only GEOMETRYCOLLECTIONs.

Either way, thanks for everyone's help!

0

from @ziggy answer,

Expand the geometrycollection with ST_CollectionExtract to query each poly inside the collection. 3 as the second argument is for Polygon, 1 and 2 are for points and lines.

ST_Intersects( ST_CollectionExtract( z.the_geom, 3 ), main.mypoint)

======

Another solution is to unfold all geometries found with ST_Dump. This will include all geometry types, not just polygons.

ST_Intersects( (ST_Dump(z.the_geom ).geom), main.mypoint)

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.