title: Python Shapely Cant Find Points Within Polygon When Looping Polygons
body: <p>I have a <strong>csv file</strong> of store locations with lat and lon. I also have a <strong>geojson</strong> file with polygon features of US Census tracts. I want to use <strong>Python</strong> to see in which polygons each location exists.</p>
  
  <p>I am using <strong>Shapely</strong> Python library's contains() to loop through the csv file of store locations, grab lat and lon coordinates, and then check if the coordinates exist within one of the polygons from a geojson file.</p>
  
  <p>The code I have works successfully if I loop through each location/coordinate first, then loop through each polygon, using contains() or within() to check if a polygon contains the point. </p>
  
  <p>But this is too slow, so I want to inverse the process and loop through the polygons/geojson file first, and check the hundreds of thousands of polygons against the thousands of coordinates, instead of the other way around.</p>
  
  <p>However, when I simply switch the loops to check the polygons against the coordinates, contains() does not find any matches, even though it finds correct matches when I loop through coordinates first and check them against polygons. It's the same code, except reversed.</p>
  
  <p>Code (contains brings up false):</p>
  
  <pre><code>with open(mapFile) as f:
      data = json.load(f)
  
      #Loop through features in Geojson file
      for feature in data['features']:
  
          polygon = shape(feature["geometry"])
  
          for row in locationsFile:
              #get point coordinates from file
              pointLat = float(row[13])
              pointLon = float(row[14])
  
              point = Point(pointLon, pointLat)
  
              print(polygon.contains(point))   
              if polygon.contains(point):
  
                  #Grab data
                  newx = feature["properties"]["x"]
                  newy = feature["properties"]["y"]
                  newz = feature["properties"]["z"]
                  #append data
                  row.append(newx)
                  row.append(newy)
                  row.append(newz)
                  #update file
                  newFile.writerow(row)
                  break
  </code></pre>
  
  <p>This code produces accurate results:</p>
  
  <pre><code>with open(mapFile) as f:
      data = json.load(f)
  
      #Loop through coordinates in CSV file
      for row in locationsFile:
          #get point coordinates from file
          pointLat = float(row[13])
          pointLon = float(row[14])
  
          point = Point(pointLon, pointLat)
  
          #Loop through features in Geojson file
          for feature in data['features']:
  
              polygon = shape(feature["geometry"])
  
              if polygon.contains(point):
  
                  #Grab data
                  newx = feature["properties"]["x"]
                  newy = feature["properties"]["y"]
                  newz = feature["properties"]["z"]
                  #append data
                  row.append(newx)
                  row.append(newy)
                  row.append(newz)
                  #update file
                  newFile.writerow(row)
                  break
  </code></pre>