Create Latitude/Longitude Arrays for GeoTIFF Image

QUESTION: Can you show me how to create latitude and longitude arrays for each pixel in a GeoTiff image?

ANSWER: This is quite easy. Details of how to read a GeoTIFF image file and how to translate the geotags built into the file as metadata can be found in the article Overlaying Geographic Data on GeoTIFF Images.

So you can play along if you like, I have made a GeoTIFF image available to you on my web page. Put the GeoTIFF image to the same directory in which you will run the IDL code in this article.

To open the GeoTIFF file and read it into IDL, type these commands.

   filename = 'AF03sep15b.n16-VIg.tif'
   image = Read_Tiff(filename, GEOTIFF=geotag)
   image = Reverse(image, 2)
   s = Size(image, /DIMENSIONS)

As described in the previous article, you can use the tie points and the scale factors in the GeoTIFF file to calculate the corners of the image. Remember, though, that these corner points are on the outside edge of the image. More typically, when you are create latitude and longitude arrays for each image pixel, you want to value at the center of each pixel. Thus, you will have to move the corner point half a pixel in U and V to make it represent the center of the corner pixel. Here is code that moves the origin from the upper left corner to the bottom left, to correspond with IDL's image origin.

   xscale  = geotag.ModelPixelScaleTag[0]
   yscale  = geotag.ModelPixelScaleTag[1]
   tp      = geotag.ModelTiePointTag[3:4]
   xOrigin = tp[0]
   yOrigin = tp[1] - (yscale * s[1])

Now, we have to add one-half pixel in both the U and V directions to center this point in the origin pixel.

   xOrigin = xOrigin + (xscale * 0.5)
   yOrigin = yOrigin + (yscale * 0.5)

Next, we create a vector of U and V values.

    uvec = Findgen(s[0]) * xscale + xOrigin
    vvec = Findgen(s[1]) * yscale + yOrigin

Next, we want to do some dimensional juggling tricks to create 2D arrays of these map coordinate values for later transformation to longitude and latitude values. (A previous version of this article left this step out and created erroneous values. Thanks to Allard de Wit for pointing this out to me.)

   uarray = Rebin(uvec, s[0], s[1])
   varray = Rebin(Reform(vvec, 1, s[1]), s[0], s[1])

Next, we set up the map projection structure. As shown in the previous article, this image is in a Albers Equal Area Conic projection, so we set the map projection structure up like this. When you do this yourself, of course, you will set up the map projection structure according to whatever map projection information is present in the GeoTiff file. The information you need to do so is contained in geotag structure returned from the file. Since this is always a pain to figure out, I just get a map coordinate object, from which I can obtain the map structure, by using the Coyote Library routine cgGeoMap. The cgGeoMap program can parse most of the common map projections.

    mapStruct = MAP_PROJ_INIT('Albers Equal Area Conic', $
       DATUM=8, $ ; WGS84
       CENTER_LAT=geotag.PROJNATORIGINLATGEOKEY, $
       CENTER_LON=geotag.PROJNATORIGINLONGGEOKEY, $
       STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
       STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)

Or, using cgGeoMap , like this.

   mapCoord = cgGeoMap(filename)
   mapStruct = mapCoord -> GetMapStruct()

And then we inverse transform the U and V arrays from UV space to latitude/longitude space.

    lonlat = MAP_PROJ_INVERSE(uarray, varray, MAP_STRUCTURE=mapStruct)    
    lon = Reform(lonlat[0,*], s[0], s[1])
    lat = Reform(lonlat[1,*], s[0], s[1]) 

Now we have latitude and longitude arrays for each pixel in the image.