1

I have a shapefile with a bunch of wildfire perimeter data, and a dataframe with lat and long columns about the lat and long of fire origin points. The information is not connected except by date, so for each fire perimeter, I am trying to check which origin points occurred on that day and fall inside the shape. To do so I open the csv files as dataframes o and p, and then I open the shapefile with

driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open(r'C:\Users\Main\Desktop\locs\here\firep18_1.shp')

# Get Projection from layer
layer = shp.GetLayer()
spatialRef = layer.GetSpatialRef()
print (spatialRef)

# Get Shapefile Fields and Types
layerDefinition = layer.GetLayerDefn()
print(layerDefinition)
inFeature = layer.GetNextFeature()

Then I run the following code.

while inFeature:
    shape = inFeature.GetGeometryRef()
    print(shape)
    break
    """ get the date attribute for the input feature """
    date = inFeature.GetField('ALARM_DATE')
    date = pd.to_datetime(date)    

    if date:
        if int(inFeature.GetField('CAUSE')) == 14:
            incidents.append(inFeature.GetField('INC_NUM'))
            temp = o[o['ALARM_DATE'] == date] 

            lat, long = list(temp.LAT83), list(temp.LON83)
            for i in range(len(lat)):
                x = lat[i]
                y = long[i]
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(x,y)
                print(x, y, point, sep='\t')
                if shape.Contains(point):
                    print("something")
                    break

    # destroy the input feature and get a new one
    inFeature = None
    inFeature = layer.GetNextFeature()

However, the point is not being created correctly on lines

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x,y)

I have tried creating a small contrived example here

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-121.4347, 39.8134)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddPoint(-122, 40)
poly.AddPoint(-120,40)
poly.AddPoint(-120,39)
poly.AddPoint(-122,39)

print(type(poly))
print(poly.IsEmpty())
print(poly.Contains(point))

which also does not work, even though I am sure the point should be within the polygon. Could anyone help me understand what AddPoint takes as arguments, because I think giving it lat long in the form (-121.4347, 39.8134) is incorrect, but I could not find useful examples or documentation.

Edit: When I run shape = inFeature.GetGeometryRef(); print(shape) I get POLYGON ((138036.905699998 -402646.3633,138086.986 -402651.0541...138036.905699998 -402646.3633)) But I am not sure how they get these values. At first I thought it was just lat and long * 1000, but the multiplied value is different for long and lat.

  • 1
    because I think giving it lat long in the form (-121.4347, 39.8134) is incorrect, well it can't be lat/long as lat has values +/- 90, so looks like you have long/lat – nmtoken Nov 8 '19 at 15:20
  • What projection is the shapefile in? – nmtoken Nov 8 '19 at 15:20
  • @nmtoken I have added an edit that shows the format of one of the shapes in the shapefiles, you are right that it is long/lat, it was a typo by me. As you can see from the shape of the shapefile object, the points are stored as long/lat. Also I apologize but I do not know what you mean by projection, this is my first time working with GIS data. – Connor Hale Nov 8 '19 at 16:10
  • This would be very easy to do in postgis using sql. Create two tables, then query using st_within(point,poly) and date1 = date2 – vinh Nov 8 '19 at 16:35
  • units look like metres (or possibly feet), not degrees. What does print(spatialRef) give you, or what is the content of the .prj file that comes with the shapefile – nmtoken Nov 8 '19 at 19:12
1

From the Python GDAL/OGR Cookbook: Create a Polygon (and the correct formulation is Point(x,y) -> Point(longitude, latitude) )

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-122, 40)
ring.AddPoint(-120, 40)
ring.AddPoint(-120, 39)
ring.AddPoint(-122, 39)
ring.AddPoint(-122, 40) # close the LinearRing - Polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
print(poly.Contains(point))
True

But with the 4 points only

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-122, 40)
ring.AddPoint(-120, 40)
ring.AddPoint(-120, 39)
ring.AddPoint(-122, 39)
poly = ogr.Geometry(ogr.wkbPolygon)
print(poly.Contains(point))
ERROR 1: IllegalArgumentException: Points of LinearRing do not form a closed linestring
False

It is easier with Shapely

from shapely.geometry import Point, Polygon
point = Point(-121.4347, 39.8134)
List_points = [Point(-122,40),Point(-120,40),Point(-120,39),Point(-122,39)]
poly = Polygon([[p.x, p.y] for p in List_points])
print(poly.contains(point))
True

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.