Thursday, January 30, 2014

Selecting and Downloading Raster Data via ftp using Python

My problem was to browse through a large number of Radarsat-2 images on a ftp-server and download the ones falling into my area of interest. Each image is available as a *.zip file and a corresponding *.tif map projected quicklook

There is two main pieces that help solving this problem:

1) You can in fact use gdalinfo to retrieve information about a raster file on ftp without having to download the file itself. The command (which I found documented here ) is, the output is here redirected to the output_tempfile:

gdalinfo /vsicurl/ftp://username:password@ftp.adress.no/outfolder/filename.tif >  output_tempfile


or calling from within Python using os.system with the variables defined before:

os.system(r'gdalinfo /vsicurl/ftp://' + username + ':' + password + "@" + ftpserver + workingfolder + file + ' > ' + output_tempfile)


2) The second step is then to read the centre coordinate from the output_tempfile and if matching download the file. The ftp download is described here , let's see the whole script, which hopefully is self-explanatory with its comments:


import ftplib, os


#DEFINE FTP SERVER PARAMETERS

ftpserver = 'ftp.youradress.no'
username = 'anonymous'
password = 'guest'
workingfolder = '/Out/max/'

#temporary file storing location info from quicklook
output_tempfile = 'C:\Users\max\Desktop\\test.txt'
#destinationfolder for download
destination_folder = 'C:\Users\max\Documents\\'

#Open ftp
ftp = ftplib.FTP(ftpserver, username , password)
ftp.cwd(workingfolder)

#Create a file list of all quicklooks
filelist = []

try:
    filelist = ftp.nlst('*.tif')
except ftplib.error_perm, resp:
    if str(resp) == "550 No files found":
        print "No files in this directory"
    else:
        raise

for file in filelist:
    
    #get quicklook file info and save in temporary file
    os.system(r'gdalinfo /vsicurl/ftp://' + username + ':' + password + "@" + ftpserver + workingfolder + file + ' > ' + output_tempfile)
    
    #Read from temporary file the centre coordinates of the quicklook    
    search = open(output_tempfile)

    for line in search:
        if "Center" in line:
            longitude_quicklook =  float(line[15:25])
        
            latitude_quicklook =  float(line[28:38])
        
        
    search.close()

    print 'Longitude ', longitude_quicklook
    print 'Latitude ', latitude_quicklook       


    #Check if scene is contained in area of interest, divided in two areas since not rectangular
    if ((-23.0 < longitude_quicklook < 67.0) and (71.0 < latitude_quicklook < 90.0)):
        contained = True
        print 'true'
    elif ((-33.0 < longitude_quicklook < 14.0) and (66.0 < latitude_quicklook < 71.0)):
        contained = True
        print 'true'
    else:
        contained = False

 
    zipfile = file[:-7] + '.zip'
    zipsavefile = destination_folder + file[:-7] + '.zip'
    
    print zipfile
    print zipsavefile
    
    #If file contained in area of interest, download the file
    if contained == True:
        print 'transferring ', zipfile
        ftp.cwd(workingfolder)
        ftp.retrbinary('RETR ' + zipfile, open(zipsavefile, 'wb').write)
        print 'transferred ', zipfile
        
    #Close ftp connection    
    ftp.quit()
        
#Remove temporary file
os.remove(r'C:\Users\max\Desktop\test.txt')

#End

Tuesday, January 21, 2014

Importance of applying SAR terrain correction

There are still many publications with SAR images which are not terrain corrected. In SAR images there is a displacement of the pixel dependent on its elevation. A few years ago, applying terrain correction was a major task, but the NEST Toolbox allows that to be done very easily by anyone (see the four post describing this here )

A very good example on the importance of applying terrain correction I had just today. The following image was geocoded simply using gdalwarp as explained in previous post.. It shows displacement for a vehicle  track visible on the SAR image in the snow compared to its true location Note that the track is displaced the same direction as the SAR layover of the mountain:



Applying Range Doppler terrain correction using NEST and a not very perfect DEM, the track matches much better:



(Subset of Antarctis Radarsat-2 image, provided by NSC/KSAT under the Norwegian-Canadian Radarsat agreement 2013)

Tuesday, January 14, 2014

Raster Calculations with numpy.where

The following may seem obvious to experienced users, but I did not find much documentation on the web on this:

With satellite images, one typically wants to do raster calculations. The easy way is to loop through the raster pixel by pixel, but this takes a long time with huge rasters. In this case, I want to calculate the precentage of sea ice. The daily percentage of the input "iceraster" is cumulatively added to the outputraster:


    #Loop through all files to do calculation
    for infile in filelist:
  
        (infilepath, infilename) = os.path.split(infile)
        print 'Processing ', infilename
  
        #open the IceChart
        icechart = gdal.Open(infile, gdalconst.GA_ReadOnly)
        if infile is None:
            print 'Could not open ', infilename
            return
  
        #get image size
        rows = icechart.RasterYSize
        cols = icechart.RasterXSize
          
        #get the bands 
        outband = outraster.GetRasterBand(1)
      
        #Read input raster into array
        iceraster = icechart.ReadAsArray()
  
        #Process the image by looping through the pixels
        for i in range(rows):
            for j in range(cols):
                if iceraster[i,j] == 0:
                    outarray[i,j] = outarray[i,j] + (  0.0 / NumberOfDays) 
                elif iceraster[i,j] == 5:
                    outarray[i,j] = outarray[i,j] + (  5.0 / NumberOfDays) 
                elif iceraster[i,j] == 25:
                    outarray[i,j] = outarray[i,j] + ( 25.0 / NumberOfDays)
                elif iceraster[i,j] == 55:
                    outarray[i,j] = outarray[i,j] + ( 55.0 / NumberOfDays)
                elif iceraster[i,j] == 80:
                    outarray[i,j] = outarray[i,j] + ( 80.0 / NumberOfDays)
                elif iceraster[i,j] == 95:
                    outarray[i,j] = outarray[i,j] + ( 95.0 / NumberOfDays)
                elif iceraster[i,j] == 100:
                    outarray[i,j] = outarray[i,j] + (100.0 / NumberOfDays)
          
              
                if iceraster[i,j] == 999:
                    outarray[i,j] = 999.0


The better way is to do calculations on the whole array/raster to use numpy.where. It is in a way similar to the band math-tool in some software packages. The syntax is:


numpy.where ( [condition], [if TRUE do this], [if FALSE do this] )

Using this, the if-loop above can be rewritten like this (and does the calculation VERY much faster):

        
outarray = numpy.where((iceraster==0.0),(outarray + ( 0.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==5.0),(outarray + ( 5.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==25.0),(outarray + (25.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==55.0),(outarray + (55.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==80.0),(outarray + (80.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==95.0),(outarray + (95.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==100.0),(outarray + (100.0/NumberOfDays)), outarray)
outarray = numpy.where((iceraster==999.0), 999.0 , outarray)


The code above checks again the iceraster value and does the calculation, otherwise the output array is laft as it is.

The output is monthly ice percentage as seen here: