Monday, March 4, 2013

Installing Python with gdal on a Windows computer

Getting a new computer gave me the chance to try again installing Python with gdal. There may be various ways, but this works for me. It can be a bit tricky to get all working, so these resources helped me a lot:

My Computer: Windows 7, 64-bit

1) Download and install Python 2.7 (the 3er versions do not support all libraries yet (?) and my code is all 2.7)

2) For installing gdal I follow these instructions. This installs both the command line gdal utilities and the python bindings to use gdal within python. I need the command line utilities also to call them with "os.system" from within python.

[A day later I realize -- the bindings from the above mentioned webpage link the command line installation to python. For some reason in my new computer I don't get this to work. Solving that quickly and being lazy: Install the Python gdal bindings from the link at point 3 below. The difference is that these bindings contain the whole gdal. BUT if you only install these, you won't have the command line utilities available]

3) At this great page, "all" libraries are available in different version, especially Windows 64-bit:
  I use it for installing additional libraries like Numpy, Pil, etc

4) I use Spyder for writing code (looks a bit like Matlab and works really great!). First install the Python GUI library "PyQt 4" from http://www.riverbankcomputing.com/software/pyqt/download and then Spyder, available here https://code.google.com/p/spyderlib/ or at http://www.lfd.uci.edu/~gohlke/pythonlibs/

[today, 1/3/2013, the 4.10 PyQt version had an issue with Spyder but that should be resolved soon in the download]






Wednesday, January 23, 2013

Geocoding Radarsat-2 images with gdalwarp

Trying to use gdalwarp to geocode a Radarsat image for creating map projected quick looks, I used the imagery_HH.tif file only, but the result was not good:


After a long search I realized that not the whole geolocation grid was accessed using imagery_HH.tif only.

To my surprise I realized that gdalwarp actually can process the Radarsat product.xml file. This file can be given to gdalwarp as inputfile, and gdalwarp will process all available bands and use the whole geolocation grid follwoing with the data.

The whole code loops through the zip file, unzips them and geocodes the image, then removing the unzipped files, leaving the original zip-file:

# -*- coding: utf-8 -*-
"""
Created on Wed Jan 16 15:43:07 2013
Takes the Radarsat-2 zipfile and creates a map projected quicklook from the imagery_HH file
@author: max """ import zipfile, glob, os, shutil #List of zip-files to be extracted #filelist = glob.glob(r'C:\Users\max\Documents\trash\*.zip') filelist = glob.glob(r'Z:\Radarsat\Sathav\2013\01_Januar\*.zip')
#Loop through zipfiles, extract and create quicklook for radarsatfile in filelist:
    #Split names and extensions     (infilepath, infilename) = os.path.split(radarsatfile)     (infileshortname, extension) = os.path.splitext(infilename)
    #Open zipfile     zfile = zipfile.ZipFile(radarsatfile, 'r')     print     print "Decompressing image for " + infilename + " on " + infilepath
    #Extract imagery file from zipfile     imagefile = infileshortname + '/imagery_HH.tif'     zfile.extractall(infilepath)
    #Define names
    #gdalsourcefile = infilepath + '\\' + infileshortname + '\\imagery_HH.tif'     gdalsourcefile = infilepath + '\\' + infileshortname + '\\product.xml'     outputfilename = infilepath + '\\' + infileshortname + '\\' + infileshortname + '_EPSG3575.tif'

    #Call gdalwarp     print     print "map projecting file"     print     os.system('gdalwarp -tps  -t_srs \"+proj=laea +lat_0=90 +lon_0=10 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs\" ' + gdalsourcefile + ' ' + outputfilename )

    #Call gdaltranslate     print     print "downsampling file"     print     os.system('gdal_translate -ot byte -outsize 8% 8% -scale 0 50000 0 255 ' + outputfilename + ' ' + browseimage )
    #Remove folder where extracted and temporary files are stored     shutil.rmtree(infilepath + '\\' + infileshortname )
    #Close zipfile     zfile.close()

print   print 'Done creating quicklooks for ', infilepath
#end 


Now the result is ok as wanted: