[gdal-dev] Rotating a GeoTransform from a KML Superoverlay

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Rotating a GeoTransform from a KML Superoverlay

Chris Moulton
Hello GDAL Dev community!

First time question asker here, and it doesn't seem like this specific task has been asked about before.

Problem: I have a bunch of KML superoverlays which I am tasked with converting them all to geotiff. I am running into an issue with the rotation values not being 0.0, which results in improperly georeferenced rasters. I've attempted a few different solutions, which I'll outline below.

Here is a sample KML bounding box (exact coords omitted for anonymity):

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<GroundOverlay>
     <name>KML Name</name>
     <Icon>
          <href>image_path.jpg</href>
          <viewBoundScale>0.75</viewBoundScale>
     </Icon>
     <LatLonBox>
          <north>40.5</north>
          <south>40.4</south>
          <east>-105.2</east>
          <west>-105.1</west>
          <rotation>58.0</rotation>
     </LatLonBox>
</GroundOverlay>

The code snippet below is my current iteration, the only inputs are IN_KML_PATH, OUT_TIFF_NAME, and then the ROTATION value in degrees read directly from the KML LatLonBox.

Here it is:

'''
#Dependencies
import gdal
from affine import Affine

#Driver for writing geotiff
tiff_driver = gdal.GetDriverByName('GTiff')

#Setting up data source, getting GeoTransform and other relevant info
in_ds = gdal.OpenEx(IN_KML_PATH)
in_geotransform = in_ds.GetGeoTransform()
in_projection = in_ds.GetProjection()
primary_band = in_ds.GetRasterBand(1)

#Creating an affine geotransform and rotating it 58.0 degrees anticlockwise (KML specs)
#This is obviously where I'm going wrong, but don't quite know what I'm missing!
af_geotransform = Affine.from_gdal(*in_geotransform)
af_geotransform *= Affine.rotation(-58.0)

#Creating out tiff 
out_ds = tiff_driver.Create(OUT_TIFF_NAME, 
     primary_band.XSize, 
     primary_band.YSize, 
     in_ds.RasterCount, 
     primary_band.DataType)

#Setting projection and new geotransform
out_ds.SetProjection(in_projection)
out_ds.SetGeoTransform(af_geotransform.to_gdal())


#Writing in bands to the output image
for i in range(1, in_ds.RasterCount + 1):
     in_band = in_ds.GetRasterBand(i)
     in_data = in_band.ReadAsArray()
     out_band = out_ds.GetRasterBand(i)
     out_band.WriteArray(in_data)

#Fin
out_ds.FlushCache()
del out_ds
'''

I've tried a few different methods for rotating the in_geotransform, including:directly changing the 'rotation' values from default 0.0 to the rotation value in the KML. This does not result in anything resembling what I need. I assume I'm missing a translation step in the affine transformation, but haven't gotten too far there. 

I know this was a lengthy question, I tried to be as detailed as possible about what I've attempted so far. I don't think there is anything going wrong per say, just reaching out to the community to see if I could get any direction here! 


Cheers,

Chris Moulton

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev