[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
cmoultonatx at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200227/5b73e76f/attachment.html>