Look at the result of evaluating
RegionUnion[p1, p2]
It is of the form
RegionUnion[Polygon[{GeoPosition[{40.6901, -73.9}], GeoPosition[{40.69, -73.8979}],...]
Polygon
requires a list of coordinates, the replacement rule /. GeoPosition -> Identity
removes the GeoPosition
GeoPosition[{40.6901, -73.9}] /. GeoPosition -> Identity
(* {40.6901, -73.9} *)
So the final result is the union of the polygons. Another way to do it
RegionUnion[p1, p2] /. GeoPosition[x_] :> x