Polygon intersection error in Shapely: "shapely.geos.TopologicalError: The operation 'GEOSIntersection_r' produced a null geometry"

As previously pointed out, p1 is not valid. On plotting it, I noticed a little 'bowtie' at the lower right. I assume you don't need this in your polygon; if not, you can try Shapely's buffer(0) trick (documented in the Shapely Manual) to fix that:

In [382]: p1.is_valid
Out[382]: False

In [383]: p1 = p1.buffer(0)

In [384]: p1.is_valid
Out[384]: True

buffer(0) has the following effect:

Before:

After:

And you can now do this:

print p1.intersection(p2)
POLYGON ((34.9396324323625151 -88.1614025927056559, 34.8937999999999988 -88.1020999999999930, 34.7998910000000024 -88.1137513649788247, 34.7998910000000024 -87.2202139999999986, 34.9994660069532983 -87.2202139999999986, 35.0041000000000011 -88.1954999999999956, 34.9396324323625151 -88.1614025927056559))

Note that I have had some instances (with areas that looked more like "birds' nests" than simple bowties) where this didn't work; check to make sure you get back a single Polygon object and not a MultiPolygon one.

You are getting this exception because p1 is not a valid polygon.

>>> p1.is_valid
False
>>> p2.is_valid
True

The documentation says that:

A valid Polygon may not possess any overlapping exterior or interior rings.

Keep in mind that since the first and last point of your polygons are different shapely is going to append the first point to the end of the list.

>>> list(p1.exterior.coords)
[(35.004100000000001, -88.195499999999996), (34.991799999999998, -85.606800000000007), (32.840400000000002, -85.175600000000003), (32.259300000000003, -84.892700000000005), (32.153500000000001, -85.034199999999998), (31.794699999999999, -85.135800000000003), (31.52, -85.043800000000005), (31.3384, -85.083600000000004), (31.209299999999999, -85.106999999999999), (31.002300000000002, -84.994399999999999), (30.9953, -87.600899999999996), (30.942299999999999, -87.592600000000004), (30.853899999999999, -87.625600000000006), (30.674499999999998, -87.407200000000003), (30.4404, -87.368799999999993), (30.1463, -87.524000000000001), (30.154599999999999, -88.386399999999995), (31.893899999999999, -88.474299999999999), (34.893799999999999, -88.102099999999993), (34.947899999999997, -88.1721), (34.910699999999999, -88.146100000000004), (35.004100000000001, -88.195499999999996)]

The exterior of your polygon is a linear ring, it also appears to be invalid:

>>> p1.exterior.type
'LinearRing'
>>> p1.exterior.is_valid
False

You can also see that if you were to turn the exterior of the polygon into a linestring it would be complex:

>>> l1 = LineString(p1.exterior.coords)
>>> l1.is_simple
False

Somehow the exterior of your polygon crosses or touches itself.

Digging a bit more into the data, we can visualize it on a map:

>>> import cgpolyencode
>>> encoder = cgpolyencode.GPolyEncoder()
>>> encoder.encode((y, x) for x, y in p1.exterior.coords)
{'points': 'svstEzthyOzkAkrxNfecL_fsAznpBcgv@ftSjsZnaeA~yRzst@_~P~mb@vwFzeXfqCvlg@w~Tvj@ra|NfjI{r@ngPfmEf`b@_ti@bvl@_oFbmx@~h]{r@~lgDsurIjdPk|hQgugAaqIntLlgFoaDwfQvsH', 'numLevels': 18, 'zoomFactor': 2, 'levels': 'PPLMKMKGKPNIKLMNPLLKJP'}

If you plug that into Google's Polyline Encoder you can see that it is the state of Alabama. If you zoom into the top left portion of Alabama you can see that two of the lines cross each other. To fix this you either need to remove the last point or swap the last point with the second to the last point.