0

I'm trying to locate the observations that lie within a collection of Polygons but my current implementation takes forever, and I suspect that it is because the 'within' method needs to check all polygons.

My current setup looks the following:

# A snip of the data, which is stored in a GeoDataFrame named
# reducedGeometryFrame in the following code.
    geometry    Centroids
0   POINT Z (-127.688834 70.55800000000001 0)   POINT (-127.688834 70.55800000000001)
1   POINT Z (-133.714 70.0585 0)    POINT (-133.714 70.0585)
2   POINT Z (-133.717167 70.0585 0) POINT (-133.717167 70.0585)
3   POINT Z (-158.410434 72.7996 0) POINT (-158.410434 72.7996)
4   POINT Z (-158.702167 72.615167 0)   POINT (-158.702167 72.615167)
5   POINT Z (-158.412667 72.800167 0)   POINT (-158.412667 72.800167)
6   POINT Z (-139.020667 70.432334 0)   POINT (-139.020667 70.432334)
7   POINT Z (-135.011 70.86966700000001 0)  POINT (-135.011 70.86966700000001)
8   POINT Z (-135.019167 70.869 0)  POINT (-135.019167 70.869)
9   POINT Z (-161.500184 71.600117 0)   POINT (-161.500184 71.600117)
10  POINT Z (-161.52755 71.599767 0)    POINT (-161.52755 71.599767)
11  POINT Z (-160.166667 76 0)  POINT (-160.166667 76)
12  POINT Z (-150 75.000167 0)  POINT (-150 75.000167)
13  POINT Z (-139.9825 73.994 0)    POINT (-139.9825 73.994)
14  POINT Z (-150 77.99983400000001 0)  POINT (-150 77.99983400000001)
15  POINT Z (-168.952167 66.3265 0) POINT (-168.952167 66.3265)
16  POINT Z (-168.568 65.781334 0)  POINT (-168.568 65.781334)
17  POINT Z (-168.263 65.746 0) POINT (-168.263 65.746)
18  POINT Z (-127.660334 70.576667 0)   POINT (-127.660334 70.576667)
19  POINT Z (-126.871334 70.681167 0)   POINT (-126.871334 70.681167)
20  POINT Z (-94.984167 69.51266699999999 0)    POINT (-94.984167 69.51266699999999)

My current implementation to separate the observations:

def separatedata():

    start = time.time()
    subSeas = [navareaIV,navareaXII,hydroArc,hydroLant,hydroPac]
    subSeasObservations = [{} for i in np.arange(len(subSeas))]
    printRange = np.arange(0,reducedGeometryFrame.shape[0],500)
    for i,obs in enumerate(reducedGeometryFrame['Centroids']):
        if i in printRange:
            print('This is the %ith iteration' % i)

        p1 = Point(list(obs.coords))

        for j,subsea in enumerate(subSeas):
            if any([p1.within(poly) for poly in subsea['geometry']]):
                subSeasObservations[j][i] = 0 

    end = time.time()

    print('It took %f second to divide the data into the five subgroups.' % (end-start))

The implementation is supposed to separate the data within five categories, based on the NAVAREA Warnings Areas, found here: https://msi.nga.mil/MSISiteContent/StaticFiles/Images/navwarnings.jpg

But those areas do not appear to be available in a shapefile online, so I constructed them myself, based on a shapefile containing the world seas, found here: http://www.marineregions.org/downloads.php (IHO Sea Areas Version 3).

The five areas of interest are constructed as follows, based on the worldSea shapefile:

worldSeas = gpd.read_file('World_Seas_IHO_v3/World_Seas_IHO_v3.shp')

navareaIV = worldSeas.loc[[18,19,20,21,26,27,28,31,32,33,37,39,40,41,43,44,45,55,56,57,58,74 
,75,76,77,78,79,80,81,84,87,89,90,95,97,99]]
navareaIV = navareaIV.reset_index(drop=True)
navareaXII = worldSeas.loc[[17,22,23,51,52,53,54,72,73,91,92]]
navareaXII = navareaXII.reset_index(drop=True)
hydroLant = worldSeas.loc[[0,61,100]]
hydroLant = hydroLant.reset_index(drop=True)
hydroPac = worldSeas.loc[[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,34,35,36,38,42,46,47,48,49,50,63,64,65,66,67,68,69,70,71,88,96]]
hydroPac = hydroPac.reset_index(drop=True)
hydroArc = worldSeas.loc[[24,25,29,30,59,60,62,82,83,85,86,93,94,98]]
hydroArc = hydroArc.reset_index(drop=True)

The above subsetting is what I expect to be the reason for why the division of the observations takes forever. For that reason, I tried to come up with a way to limit the number of polygons, that was needed to check in terms of if a point is contained within it. The reduction technique I have tried was to construct a Polygon of the union of the boundaries of each of the polygons with a certain area, and then check if a point is located with that new polygon. However, this gives me a new problem, as the polygon does not seem to contain a point, that obviously lies within the polygon, using the 'within' method and the 'contain' method.

# hydroArc example of the reduction technique
seaindecies = [24,25,29,30,59,60,82,83,85,86,93,94,98]# This is the northen part of the area, as 62 contains all of the southen part.

unionList = [worldSeas['geometry'].loc[subsea].boundary for subsea in seaindecies]
hydroArcArea = unary_union(unionList)

# Creating a geoDataFrame based on the new union boundaries.
hydroArcNorth = gpd.GeoDataFrame({'geometry':hydroArcArea},geometry='geometry')

# The area can be plotted in the following way, to see that the proposed point obviously lies within the area
hydroArcNorth.plot(figsize=(10,10))
plt.xlim(-190,190)
plt.ylim(60,95)
plt.show()

# Point 
p1 = Point((0.0,89.0))
# Checking if the point is within the area
p1.within(hydroArcNorth)

This returns the error: 'GeoDataFrame' object has no attribute '_geom'. So, instead I tried to see if the 'hydroArcNorth' contained the point:

any(hydroArcNorth.geometry.contains(p1))

Which just returns false.

I seek help in terms of optimizing the separation function above, or why the 'within' method returns the error it does or perhaps a completely different way of tackling this problem.

  • Have you tried sjoin with option within? – BERA Apr 1 at 14:08
  • @BERA Do you mind clarifying your suggestion a bit? – Kristian Nielsen Apr 1 at 15:28
  • 1
    I have a better, in terms of running time, solution incorporating 'sjoin', and I'll post it for future needs as soon as I'm done. – Kristian Nielsen Apr 2 at 12:51

Your Answer

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

Browse other questions tagged or ask your own question.