Create a point from a lat/long list using GDAL

Question

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.

Answer 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