Internal: GeoTIFF to KML example
Anonym
TOPIC:
This help article shows an example of how to read a GeoTIFF file's contents into IDL, generate an image of the contents using the IMAGE function, and then use the SAVE method to write the image to a KML file. The file used for the example is the "boulder.tif" available in the "examples/data" directory in the IDL distribution.
Discussion:
To read the "boulder.tif" into IDL and output a KML of the file, the following steps were followed:
1) The READ_TIFF routine was used with the GEOTIFF and ORIENTATION keywords to read the data from "boulder.tif" into IDL.
2) The ROTATE method was used to correct the orientation of the image.There is a table in the ORIENTATION keyword section of the READ_TIFF help page.
https://www.exelisvis.com/docs/READ_TIFF.html
3) The contents of the GEOTIFF keyword were used to determine the meaning of the GeoTIFF tag information. Information about the meaning of the GeoTIFF tags can be found at the following web address.
http://www.remotesensing.org/geotiff/spec/geotiffhome.html
4) The GeoTIFF tags ModelTiePointTag and ModelPixelScaleTag were used to determine the range of the image.
5) The MAP function was used to generate a UTM Zone 13 map.
6) To make it easier to warp the TIFF image, the MAPINVERSE method of the MAP object was used to convert Easting and Northing values to latitude and longitude.
7) The IMAGE function was then used with the IMAGE_LOCATION, IMAGE_DIMENSIONS, and OVERPLOT keywords to warp and display the TIFF image on the map generater in step 5.
8) The LIMIT of the map was reset to the values determined in step 6.
9) The SAVE method was used to save the image as a KML file.
The code to do this is shown below (with the steps labeled):
pro geo_tiff_2_KML_example
compile_opt idl2
;***************************
; Step 1
;***************************
dataFilePath = FILEPATH('boulder.tif',$
subdir=['examples','data'])
data_variable = read_tiff(dataFilePath, GEOTIFF=GeoKeys, ORIENTATION=a)
;Get the size of the IMAGE
dims = size(data_variable,/DIMENSIONS)
;****************************
; Step 2
;****************************
;Need to adjust orientation of image. The Orientation
;of the image is 1. Therefore, I need to use the ROTATE
;function with the DIRECTION argument set to 7. See the
;ORIENTATION section of the READ_TIFF help for
;more information.
data_variable = rotate(data_variable,7)
;****************************
;Step 3
;Ouput of "help, GeoKeys"
;at the bottom of the code
;****************************
;****************************
;Step 4
;****************************
;Calculate lowest Northing value for image
low_north = geokeys.ModelTiePointTag[4]- $
(dims[1]*geokeys.MODELPIXELSCALETAG[0])
;Calculate highest Easting value for image
largest_easting = geokeys.ModelTiePointTag[3] + $
(dims[0]*geokeys.MODELPIXELSCALETAG[0])
;*****************************
;Step 5
;*****************************
;Create a MAP with the UTM Zone 13
m = MAP('UTM', zone=13)
;*****************************
;Step 6
;*****************************
;Use the MapInverse method to determine the latitude and
;longitude values for the boundries of the image.
lowerleft_latlon = m.MapInverse(geokeys.ModelTiePointTag[3],$
low_north)
topright_latlon = m.MapInverse(largest_easting, $
geokeys.ModelTiePointTag[4])
;Define the DIMENSIONS of the image and LIMIT of the map
;using the LAT/LON values from the calculation above.
lim =[lowerleft_latlon[1],lowerleft_latlon[0],$
topright_latlon[1],topright_latlon[0]]
img_dims = topright_latlon - lowerleft_latlon
;******************************
;Step 7
;******************************
i = image(data_variable,$
/overplot,$ ;Use the same coordinate system as the MAP
image_location=lowerleft_latlon,$ ;Location of the lower left pixel
GRID_UNITS=2,$ ;The units of the grid are degrees
IMAGE_DIMENSIONS=img_dims) ;
;******************************
;Step 8
;******************************
;Reset the limit of the map
m.limit = lim
;******************************
;Step 9
;******************************
;Save the map to a KML file
i.save,"boulder.kml"
end
;*********************************************************************
; Information out put form GEOKEYS. Needed to look up information
; about these keys in the following location.
; http://www.remotesensing.org/geotiff/spec/geotiffhome.html
;
; MODELPIXELSCALETAG
; DOUBLE Array[3] (Softdesk) (10,10,0)
; MODELTIEPOINTTAG
; DOUBLE Array[6, 1] (Intergrah)
; [ 0.00000000 0.00000000 0.00000000
; 468034.00 4441479.0 0.00000000]
; GTMODELTYPEGEOKEY
; INT 1 (Projection Coordinate System)
; GTRASTERTYPEGEOKEY
; INT 1 (Raster pixel is Area)
; GEOGLINEARUNITSGEOKEY
; INT 9001 (EPSG Linear Units)
; GEOGANGULARUNITSGEOKEY
; INT 9102 ( Degree)
; PROJECTEDCSTYPEGEOKEY
; INT 26713 (PCS_NAD27_UTM_zone_13N)
;********************************************************************
Reviewed by DS 4/21/2014