3

I have read in a KML file of CA firezones and for a certain area (around LOMPOC, CA), it seems to include anything inside a non-fire zone doughnut hole. Most of the processing is fine, just there are a few areas surrounded by fire zone where areas outside of the zones are included as being in the fire zone areas.

Is there some way to fix the issue or use another approach?

The fire threat maps are here: https://www.cpuc.ca.gov/FireThreatMaps/

I am using the KMZ files unzipped to KML.

Code follows

import geopandas as gpd
import matplotlib.pyplot as plt
gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'

fp = "CPUC_Fire-Threat_Map.kml"
polys = gpd.read_file(fp, driver='KML')

Tier2 = polys.loc[polys['Name']=='Tier 2']
Tier2.reset_index(drop=True, inplace=True)

#test points - visually compared on map
POINT_DOUGHNUT = Point(-120.441806,34.652998)
POINT_INSIDE = Point(-120.142221,34.66407)
POINT_OUTSIDE = Point(-119.471006,36.414549)

POINT_INSIDE.within(Tier2.loc[0, 'geometry'])
> True  (correct)  

POINT_OUTSIDE.within(Tier2.loc[0, 'geometry'])
> False (correct)

POINT_DOUGHNUT.within(Tier2.loc[0, 'geometry'])
> True  (wrong -- should be false)

To add the area that seems to be the problem:

Lompoc CA area with Tier 2 firezones

Improve this question
8
  • I'm not sure the terminology but I've take a KMZ file provided by the California Public Utilities Commission, unzipped it and used the KML (document?) and read it using geopandas as above. I'm just using the Tier2 polygon from that file with the three test points. Hope this helps clarify.
    – berlin
    Nov 21, 2019 at 3:28
  • You can also go to ia.cpuc.ca.gov/firemap/# which is a rendering of the fire zones. If you search for Lompoc, CA you'll see the area that I'm concerned about. Anything in the clear Lompoc area is being found as "within" the Tier 2 firezone even thought it is a doughnut hole in the firezone.
    – berlin
    Nov 21, 2019 at 3:29
  • 1
    Are you sure the polygon is actually considered a hole or is it just an overlapping polygon. In order for the polygon to be identified as a hole it's coordinate vertice definition has to be the opposite direction of the polygon it sits in. In orther words if the start coordinate in the definition for the outer polygons moves clockwise then the internal polygon must go counterclockwise. No sure how to check this other than to try to convert it to a different fomrat WKT or shapefile and see how it looks
    – jport
    Nov 21, 2019 at 14:23
  • No, I have no idea how this was built. There are ESRI shapefiles provided as well at the link above but I'm not sure how to check this issue nor how to remedy it if it is a problem. Right now, I'll probably end up visually inspect the layouts of over 400K points in order to ensure we have our mapping correct or close. Its for a low-income electric resiliency program and it has to be spot on.
    – berlin
    Nov 21, 2019 at 14:45
  • 1
    Can you screen shot the area where the error is? I downloaded the shapefile and want to check if I can see the problem in the shapefile version.
    – jport
    Nov 21, 2019 at 14:54

2 Answers 2

Reset to default
2

It looks like there's an issue with the fact that this is a multipolygon (multiple polygons defined within a single polygon definition). In fact the entire shapefile or KMZ is defined as a single polygon with many internal rings. You might have to explode this multipolygon into separate polygons which I'm not 100% sure how to do in code (arcgis or qgis) will do this for you.

Here's a link to someone having a similar issue that you can follow up with

Check if a point falls within a multipolygon with Python

Also here's how to deal with the shapefile in ArcGIS or QGIS

Saving of multi-polygon shapefile into a individual polygon shapefiles

If you end up with strict python code on how to explode this shapefile post it back as it would be great to see for all those interested.

Improve this answer
1

The end result was to go into QGIS as suggested above and then save the multipolygon into individual shapefiles. One key to making it work was editing the attribute table to give each individual polygon shapefile its own unique FID.

Improve this answer

Your Answer

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

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