0

I'm getting an error about an operation on mixed geometries coming from a use of ST_Clip(raster, geometry). All the tables involved are uniformly SRID 4326 but maybe it's coming from an empty raster (a clip of a clip)? Any ideas for what is causing this or how to figure it out would be very welcome!

Here's the query:

select "name1", "name2",
       sum("areal_overlap_wgt" * "tract_pop") as "areal_overlap_pop",
       sum("dasymmetric_overlap_wgt" * "tract_pop") as "dasymmetric_overlap_pop",
       sum("dasymmetric_wgt" * "tract_pop") as "dasymmetric_pop"
from (
    select "name1", "name2", "tract_pop",
           case
            when "tract_area" > 0
            then "tract_area_in_both" / "tract_area"
            else 0
           end as "areal_overlap_wgt",
           case
            when "dev_in_tract" > 0
            then "dev_in_tract_s1_s2" / "dev_in_tract"
            else 0
           end as "dasymmetric_overlap_wgt",
           case
            when "dev_in_tract" > 0
            then "dev_in_tract_s1" / "dev_in_tract"
            else 0
           end as "dasymmetric_wgt"
    from (
        select "name1", "name2", "tract_pop",
               ST_AREA("tract_geom" :: geography) as "tract_area",
               ST_AREA(ST_INTERSECTION("tract_in_s1", "geom2") :: geography) as "tract_area_in_both",
               coalesce(ST_VALUECOUNT(ST_UNION("rast_in_tract"),1,true,1), 0)   as "dev_in_tract",
               coalesce(ST_VALUECOUNT(ST_UNION(ST_CLIP("rast_in_tract", "geom1", true)),1,true,1), 0)   as "dev_in_tract_s1",
               coalesce(ST_VALUECOUNT(ST_UNION(ST_CLIP(ST_CLIP("rast_in_tract", "geom1", true), "geom2", true)),1,true,1), 0)   as "dev_in_tract_s1_s2"
        from (
            select s1."name" as "name1",
                   s1."geometry" as "geom1",
                   s2."name" as "name2",
                   s2."geometry" as "geom2",
                   t."geom" as "tract_geom",
                   t."AQNFE001" as "tract_pop",
                   ST_INTERSECTION(s1."geometry", t."geom") as "tract_in_s1",
                   ST_INTERSECTION(ST_INTERSECTION(s1."geometry", s2."geometry"), t."geom") as "tract_in_both",
                   ST_CLIP(r."rast", t."geom", true) as "rast_in_tract"
            from "co_sldu" s1
            inner join "co_cd" s2 on s1."geometry" && s2."geometry"
            inner join "tracts2022_acs2017_2022" t on s1."geometry" && t."geom"
            inner join "nlcd_us_dev" r on t."geom" && r."rast"
        )
        where (not ST_ISEMPTY("rast_in_tract"))
        group by "name1", "name2", "geom2", "tract_pop", "tract_geom", "tract_in_s1"

    )
)
group by "name1", "name2"

and the error (from python which is sending that query):

Traceback (most recent call last):
  File "/Users/adam/BlueRipple/GeoData/code/areal.py", line 329, in <module>
    cur.execute(overlap_sql)
psycopg2.errors.InternalError_: contains: Operation on mixed SRID geometries (MultiPolygon, 4326) != (Point, 0)
CONTEXT:  PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 4 at IF
1
  • Try setting the srid on your inner clip ST_SETSRID(ST_CLIP(r."rast", t."geom", true),4326) as "rast_in_tract"
    – jbalk
    Commented Feb 26 at 17:14

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.