from osgeo import ogr
from osgeo import osr

def h(name):
  print "<Placemark><name>%s</name>" % (name,);

def f():
  print "</Placemark>";

print """
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>""";

gj_oslo = """{"type":"Point","coordinates":[59.95,10.75]}"""
gj_paris = """{"type":"Point","coordinates":[48.8567,2.3508]}"""

p_oslo = ogr.CreateGeometryFromJson(gj_oslo)
p_paris = ogr.CreateGeometryFromJson(gj_paris)

h("point oslo")
print p_oslo.ExportToKML()
f()
h("point paris")
print p_paris.ExportToKML()
f()

# http://geojson.org/geojson-spec.html#coordinate-reference-system-objects
# EPSG:4326
sr_geojson = osr.SpatialReference()
sr_geojson.ImportFromEPSG(4326)
# sr_geojson = p_oslo.GetSpatialReference()

# Oslo: 32V (31N: EPSG:32632)
sr_oslo = osr.SpatialReference()
sr_oslo.ImportFromEPSG(32632)

# Paris: 31U (31N: EPSG:32631)
sr_paris = osr.SpatialReference()
sr_paris.ImportFromEPSG(32631)

tr_geo_oslo = osr.CoordinateTransformation(sr_geojson, sr_oslo)
tr_geo_paris = osr.CoordinateTransformation(sr_geojson, sr_paris)

tr_oslo_geo = osr.CoordinateTransformation(sr_oslo, sr_geojson)
tr_paris_geo = osr.CoordinateTransformation(sr_paris, sr_geojson)

# If i understand correctly, UTM has a unit of 1 metre
p_oslo.Transform(tr_geo_oslo)
p_paris.Transform(tr_geo_paris)

# Buffer requires the radius to be in the same unit as the coordinates
poly_oslo = p_oslo.Buffer(600 * 1000) # 600 km
poly_paris = p_paris.Buffer(1200 * 1000) # 1200 km

# Transform back before we do the intersection
poly_oslo.Transform(tr_oslo_geo)
poly_paris.Transform(tr_paris_geo)

intersection = poly_oslo.Intersection(poly_paris)

#poly_oslo.Transform(tr_oslo_geo)
#poly_paris.Transform(tr_paris_geo)

h("poly oslo")
print poly_oslo.ExportToKML()
f()
h("poly paris")
print poly_paris.ExportToKML()
f()

intersection.Transform(tr_oslo_geo)

h("poly intersect")
# Why are you empty????
print intersection.ExportToKML()
f()

print """</Document>
</kml>"""