[GIS] Create a ground overlay kml from georeferenced raster

arcgis-10.0gdalgoogle earthkmlr

I have a world georeferenced raster (5m resolution, ~4000 x 4000 pixels) that I would like to display in google earth. I am looking for a simple and effective way to generate this kml file.

Using ArcGIS layer to KML outputs a kmz and resample the raster (I went as far as 3000 dpi, but quality is still pretty low)
Using GDAL raster2tile is slower than ArcGIS, but it also resample the image, yielding low quality output (I do not know how or if it is possible to tweak dpi or quality).
Using R maptools::GE_SpatialGrid function gave the best results so far, as it simply outputs a xml file specifying the extent of the raster. Though it does not display fully in Google Earth (no error message).

To sum it up, my preferred method would be to output a kml file, that would be entirely displayed in google earth.
I like the kml idea of simply specifying an extent, it is easy to update the raster without any more processing.

Best Answer

You need to create SuperOverlays using the GDAL tools (gdal2tiles)

Creating SuperOverlays

enter image description here

continues here:

https://developers.google.com/kml/documentation/regions#superoverlays

Here is good guide for gdal2tiles and creating KML SuperOverlays:

https://developers.google.com/kml/articles/raster

Related Solutions

[GIS] Are there sources for obtaining pre-made KML data from Google

KML is an OGC standard so there should be a lot more applications now a days which will allow you to "export to KML".

OpenLayers makes it very easy to READ and WRITE KML files using Javascript Code.

GDAL will also allow you to READ and WRITE KML files in command line. For example: create a KML file from a Spatial Database Query.

EDIT: You can also export to KML from ArcGIS 10

EDIT #2 Google does not give you the ability to download KML from it's map data. It does provide some map web services but they are not what you are thinking.

ESRI and many others do host map services for free though. For example: http://services.arcgisonline.com/ArcGIS/rest/services/Reference/World_Boundaries_and_Places/MapServer

Notice at the bottom of the above page you can "Generate KML". But what's better, this is a REST service you can query straight from code.

[GIS] Exporting 3GB ArcGIS Raster to KML without losing resolution

Save your money, buy nothing. Associate this script with a toolbox tool. Set it up with 4 parameters (Input Raster, Output Folder, Draw Order and Output Prefix). Run. Tried the SuperOverlay thing from GDAL and didn't think much of that horrible data structure. I never cleaned this up so it does processing it doesn't need to if I remember correctly, but it ran and created a working output. The kmz will even work on a Garmin (*MAP model or newer) as long as you make sure it complies with the models limitations for "Custom Maps".

#First crack at writing something to tile out a source image to
#conform to the Garmin Custom Maps KML format. For now I am not
#concerning myself with limitations of this format on various devices
#and just working on getting the output to always function in
#Google Earth
#
#To use the output kmz on most Garmin *MAP models open the kmz file in
#Google Earth. Delete tiles until there are only 100 remaining (maximum)
#Resave the resultant kmz to a new file and copy this to the "\Garmin\CustomMaps"
#folder on the GPS
#This image will display on the unit at different zoom levels based on the
#resolution of the image. If you want to use the imagery at a coarser scale
#you will have to resample the input image to a larger cell size and rerun
#this tool.
#
#Michael Stead, August 29, 2014

# Import arcpy module
import arcpy, arcgis, os, zipfile, shutil

# Get input parameters
source = arcpy.GetParameterAsText(0)
outfolder = arcpy.GetParameterAsText(1)
level = arcpy.GetParameterAsText(2)
outname = arcpy.GetParameterAsText(3)
kmzfile = outfolder + "\\" + outname + ".kmz"
workfile = "tmplayer.tif"

# Get projection of input raster... set output to wgs84
desc1 = arcpy.Describe(source)
sproj = arcpy.SpatialReference(4326)
arcpy.AddMessage(str(desc1.spatialReference.name))
arcpy.AddMessage(str(sproj.name))

#make temporary raster layer, ensuring 3 bands only
arcpy.env.workspace = outfolder
arcpy.MakeRasterLayer_management(source, workfile, "#", "#", "1;2;3")

#make sure 8 bit input as other bit depth creates empty output from Split Raster
spath = desc1.Path
if ":" in str(source):
    ddesc = arcpy.Describe(source+"\\Band_1")
if not ":" in str(source):
    ddesc = arcpy.Describe(spath+"\\"+source+"\\Band_1")
sourcedepth = ddesc.pixelType
arcpy.AddMessage(sourcedepth)
if sourcedepth <> "U8":
    arcpy.AddMessage("Convert to 8 bit unsigned")
    arcpy.CopyRaster_management(workfile, "t_raster.tif","DEFAULTS","","","","","8_BIT_UNSIGNED")
    workfile = "t_raster.tif"
else:
    arcpy.AddMessage("Already 8 bit unsigned") #but couldn't see how else to pass this as a variable
    #arcpy.CopyRaster_management ("tmplayer.tif", "t_raster.tif","DEFAULTS","","","","","8_BIT_UNSIGNED")


#check if the input is wgs84 already, if not reproject
if desc1.spatialReference.name == arcpy.SpatialReference(4326).name:
    arcpy.AddMessage("Already WGS84")
    #arcpy.CopyRaster_management ("t_raster.tif", "p_raster.tif","DEFAULTS","","","","","8_BIT_UNSIGNED")

else:
    # Reproject to WGS84
    arcpy.AddMessage("Projecting")
    arcpy.ProjectRaster_management(workfile, "p_raster.tif", sproj, "BILINEAR", "#","#", "#", "#")
    workfile = "p_raster.tif"

#Figure out size for tiles (must be less than 1024x1024)
if workfile == "tmplayer.tif":
    workfile = source
desc2 = arcpy.Describe(workfile+"\\Band_1")
if ":" in str(source):
    arcpy.AddMessage("Full Path")
    desc4 = arcpy.Describe(source+"\\Band_1")
if not ":" in str(source):
    arcpy.AddMessage("TOC Layer")
    desc4 = arcpy.Describe(spath+"\\"+source+"\\Band_1")

cellsize = desc4.MeanCellHeight
columnfactor = int((desc2.width)/1024)+1
rowfactor = int((desc2.height)/1024)+1
columnwidth = int(desc2.width/columnfactor)+1
rowheight = int(desc2.height/rowfactor)+1
arcpy.AddMessage("Height: %d" % rowheight)
arcpy.AddMessage("Width: %d" % columnwidth)
tilesize = "%d %d" % (columnwidth, rowheight)

# Create "files" folder for tiled images
arcpy.AddMessage("Creating Folder")
arcpy.CreateFolder_management(outfolder,"files")
tilefolder =  (outfolder +'\\files\\')   

# Split unprojected raster into a folder of 1024 by 1024 (or less to minimize "no data" tile area) jpg tiles
arcpy.AddMessage("Tiling Input")
arcpy.SplitRaster_management(workfile, tilefolder, outname,"SIZE_OF_TILE","JPG","CUBIC","1 1",tilesize,"0","PIXELS","#","#")

#Create empty doc.kml file
arcpy.AddMessage("Create KML")
kml_file = (outfolder +'\\doc.kml')
f = open(kml_file,"a+")
f.close()

#Add header to empty doc.kml
arcpy.AddMessage("Build KML Header")
text_file = open(kml_file, "a+")
text_file.write('<?xml version="1.0" encoding="UTF-8"?>\n')
text_file.write('<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">\n')
text_file.write('  <Document>\n')
text_file.write('    <Name>MikesKMZs</Name>\n')
text_file.close()

#Loop through folder of tiled rasters and write the extents to the doc.kml file
arcpy.AddMessage("Create KML Tile Extent")
arcpy.env.workspace = tilefolder
rasterList = arcpy.ListRasters("*", "JPG")
for raster in rasterList:
    loopath = tilefolder+'\\'+raster
    arcpy.AddMessage(loopath)
    describeimg = arcpy.Describe(loopath)
    ##### QUESTION - the west and east values might need to get swapped for eastern hemisphere, and the north and south for the southern hemisphere
    North = describeimg.extent.YMax
    South = describeimg.extent.YMin
    West = describeimg.extent.XMin
    East = describeimg.extent.XMax
    text_file = open(kml_file, "a+")
    text_file.write('    <GroundOverlay>\n')
    text_file.write('      <name>'+raster+'</name>\n')
    text_file.write('      <drawOrder>'+str(level)+'</drawOrder>\n')
    text_file.write('      <color>ffffffff</color>\n')
    text_file.write('      <Description>\n')
    text_file.write('      </Description>\n')
    text_file.write('      <Icon>\n')
    text_file.write('        <href>files/'+str(raster)+'</href>\n')
    text_file.write('      </Icon>\n')
    text_file.write('      <LatLonBox>\n')
    text_file.write('        <north>'+str(North)+'</north>\n')
    text_file.write('        <south>'+str(South)+'</south>\n')
    text_file.write('        <east>'+str(East)+'</east>\n')
    text_file.write('        <west>'+str(West)+'</west>\n')
    text_file.write('        <rotation>0</rotation>\n')    
    text_file.write('      </LatLonBox>\n')
    text_file.write('    </GroundOverlay>\n')
    text_file.close()

#Add footer to doc.kml file\
arcpy.AddMessage("Add KML Footer")
text_file = open(kml_file, "a+")
text_file.write('  </Document>\n')
text_file.write('</kml>')
text_file.close()

#Create kmz and compress kml and jpgs
zip = zipfile.ZipFile(kmzfile,'w')
zip.write(kml_file,'doc.kml',zipfile.ZIP_DEFLATED)
JPGList = arcpy.ListFiles("*.JPG")
for JPGs in JPGList:
    zip.write(tilefolder+'\\'+JPGs,'files\\'+JPGs,zipfile.ZIP_DEFLATED)

#Clean up working file debris
arcpy.env.workspace = outfolder
os.remove(outfolder+'\\doc.kml')
shutil.rmtree(tilefolder)
if arcpy.Exists("p_raster.tif"):
    arcpy.Delete_management("p_raster.tif")
if arcpy.Exists("t_raster.tif"):
    arcpy.Delete_management("t_raster.tif")    
if arcpy.Exists("tmplayer.tif"):
    arcpy.Delete_management("tmplayer.tif")

Edit - In 10.2 the line:

arcpy.SplitRaster_management(workfile, tilefolder, outname,"SIZE_OF_TILE","JPG","CUBIC","1 1",tilesize,"0","PIXELS","#","#")

Must be modified to:

arcpy.SplitRaster_management(workfile, tilefolder, outname,"SIZE_OF_TILE","JPEG","CUBIC","1 1",tilesize,"0","PIXELS","#","#")

I realized recently that you can do a find and replace for JPG to PNG in the script and it will support transparencies in Google Earth. The data I converted had 0 as "NoData" and it worked great. I haven't tested if it still works on a Garmin GPS, but their documentation says it wouldn't.

Best Answer

Related Solutions

[GIS] Are there sources for obtaining pre-made KML data from Google

[GIS] Exporting 3GB ArcGIS Raster to KML without losing resolution

Related Question