Batch Download ESA's Sentinel satellite images

The new ESA's satellite Sentinel  data is distributed free using an open data hub.
All you need to do is to register. The data hub provide tools for searching data and creating a cart of products for downloading. Since I'm very lazy, I created a Python client that would download the data once provided with the correct data link or a product list file (that can be easily downloaded from the hub online interface. The code is freely available on GitHub at:
https://github.com/rannof/SentinelDL

Next step (if time will allow it) is to add an automated data search based on parameters such as location, data type, dates etc.

Python Event Listener - multiprocessing deamon

I wanted to have a file change listener such that if new data is added to a file I could run a code.
Following this post, I've written a modification so now the listener is working in multiprocessing.

#!/usr/bin/python
# A code to stream file as it grows
# by Ran Novitsky Nof Oct 25, 2012

import os,time
from multiprocessing import Process,Pipe

# change this file name
FileName = 'stream.txt'

# adapted from http://www.valuedlessons.com/2008/04/events-in-python.html


# the Event class take care of adding, removing and executing the functions
class Event:
    def __init__(self):
        self.handlers = set()

    def handle(self, handler):
        self.handlers.add(handler)
        return self

    def unhandle(self, handler):
        try:
            self.handlers.remove(handler)
        except:
            raise ValueError("Handler is not handling this event, so cannot unhandle it.")
        return self

    def fire(self, *args, **kargs):
        for handler in self.handlers:
            handler(*args, **kargs)

    def getHandlerCount(self):
        return len(self.handlers)

    __iadd__ = handle
    __isub__ = unhandle
    __call__ = fire
    __len__  = getHandlerCount
 
# the stream class takes care of the
# function to fire in case of an event
class Stream:
  def __init__(self,fname):
    self.fname = fname
    self.fd = open(fname,'r')
    self.get_stream()
  def get_stream(self):
    # this is the actual function to fire
    print self.fd.read(),
  def restart(self):
    self.fd.seek(0)
 
# the MockFileWatcher creats the multiprocessing deamon to alert on file change
class MockFileWatcher:
    def __init__(self,source_path,sleeptime=0.1):
        self.fileChanged = Event()
        self.source_path = source_path
        self.running=False
        self.deamon=False
        self.fd = None
        self.sleep = sleeptime
        # the next two lines can be set outside the class
        self.stream = Stream(source_path)
        self.fileChanged += self.stream.get_stream # add event handler
    def watchFiles(self):
        # open the file in non-blocking read mode
        self.fd = os.open(self.source_path,os.O_RDONLY | os.O_NONBLOCK)
        # go to the end
        os.lseek(self.fd,0,2)
        self.running=True
        self.deamon=False
        while self.running:
          # try to read
          if os.read(self.fd,1)!='':
            # New data was added
            os.lseek(self.fd,0,2) # go to the end of file
            self.fileChanged() # fire the event
          # see if the parent process has a message
          if self.chiled_conn.poll():
            # execute the message
            eval(self.chiled_conn.recv())
          # sleep for a while
          time.sleep(self.sleep)
    def watchFilesDeamon(self):
        # start the deamon
        self.parent_conn,self.chiled_conn = Pipe()
        p = Process(target=self.watchFiles)
        p.start()
        self.p = p  
        self.deamon=True
        return p
    def stop(self):
        # stop deamon/process
        if self.deamon:
          self.parent_conn.send('self.stop()')
          print "Hold while stopping child process..."
          # hold depends on the sleeping time
          if self.parent_conn.recv():
            pass
            #print "chiled is stopped"
          self.p.terminate() # terminate multiprocessing
          self.deamon=False
        # stop process
        self.running=False
        self.fd=None
        # if this is the process make sure parent know we are done
        if not self.deamon:  self.chiled_conn.send(True)      
        return True

if __name__=='__main__':
  watcher = MockFileWatcher(FileName)
  p = watcher.watchFilesDeamon()
  print 'from this point the listner is running at the background'
  print 'add more code here so you can see how it works.'

Try this by changing the file name at the top and run it. then add new lines to the file.

SeisComp3 - get events data

After installing  Seiscomp3 and running it for a while, now it is time to play with the data.
As a first step I want to be able to see the events outside the SeisComp3 environment.
This will later serve me when exporting the data to a web page or reports.
The steps described here are: 1) get a list of events, 2) export the data to an xml file and 3) print the data as a bulletin.
1) Get a list of events:
scevtls -d mysql://sysop:sysop@localhost/seiscomp3
Will give a list of event codes like:

Wrong 'begin' format -> setting to None
Setting start to 1970-01-01 00:00:00
Setting end to 2012-10-24 08:17:13
gfz2012uizn
gfz2012atlf
...

where gfz will be the agency prefix and the last 4 letters will be random code.
You can also set begin and end time to get only a time window events. see scevtls documentation.

2) For this example I select one event (in csh environment):
set EVENT = gfz2012uizn
Now we can export the event data to an xml file:
scxmldump -d mysql://sysop:sysop@localhost/seiscomp3 -E $EVENT -P -M -A -o test.xml
see scxmldump for more information

3) Print the data to stdout:
scbulletin -i test.xml
Note that it is possible to skip the xmp part by using:
scbulletin -d mysql://sysop:sysop@localhost/seiscomp3 -E $EVENT
also adding -3 as a flag will output even more details.
see scbulletin for more info.

I'll update later with a python script to export the data to a Google Earth kml file.

HowTo Mosaic ASTER Global Digital Elevation Model (GDEM)

A second version of ASTER GDEM was released in Oct. 17 2011 by The Ministry of Economy, Trade and Industry of Japan (METI) and the National Aeronautics and Space Administration (NASA). The GDEM is in a 30m pixel resolution and 1x1 degree tiles. The data is distributed as zipped GeoTIFF files. When a larger than 1x1 deg. DEM is needed for Interferometry processing, you can use my GDEM.py script available from my script page. This python script will need a system with pylab and PIL python modules.

Simply:

Download the GDEM zipped files (see here on how to do that).

Unzip the zip files in the working directory.

The working directory should contain sub-directories containing the GeoTIFF files.

Run the script giving -h for help or the final DEM name as arguments.

The Final DEM is aimed for using as DEM in ROI-PAC interferometry processing.

If a -g is added as an argument, the produced DEM will be suitable for interferometric processing using Gamma software.

Create a photo collage using a single command-line

In order to create a contact-sheet or a collage of many files under linux, I prefer to use a single command-line. I create a directory with all the images and use ImageMagick montage function:

montage -font Arial-Narrow -pointsize 10 -label '%t' -resize 100x150 *.jpg -geometry +3+3 -tile 4x4 -shadow collage.jpg

RoiView - Explore InSAR data and more

RoiView was originaly designed to replace DGX/MDX software which require bureaucracy to obtain. This software was used in order to view ROI-PAC interferometry processing results but producing an image from file was not so strait forward. In time support was added for ERDAS ER Mapper header files and for user provided parameters with or without header files.

Current RoiView version (0.75) is available for download at
See also RoiView's help files

Screenshot:

HowTo Parse HTML text using python.

I needed to to get the text of a specific div elemnt in an html file. I tried to use python's standard library modules for markup processing such as htmllib etc. but I couldn't figure how to use them. I've created my own module just for getting the text from an html element:

#!/bin/env python
from htmlentitydefs import entitydefs as ent
import string
# This module enable you to extract text from a certian HTML element
#
# by Ran Novitsky Nof, 2010
# ran.nof@gmail.com
#
# example of use:
# say we want to get the text in an element of type tag (e.g. 'div','a','span' etc.)
# who has an attribute key (e.g. "id","class","href" etc.) with a value of val
# for example in order to extract the text of a div element with id of textdiv from a file htmlfile.html:
#
# <html>
# <head>
#   :
# </head>
# <body>
#   :
# <div "id"="textdiv">I will not buy this <a href="spam">record</a> it is scratched. </div>
#   :
# </body>
# </html>
#
# use:
# from htmlparser import Parser
# htmlfile='htmlfile.html'
# tag,key,val = ('div','id','textdiv')
# text=Parser(htmlfile).getText(tag,key,val)
# print(text)
#
class Element():
  def __init__(self):
    self.startTag = -1
    self.endTag = -1
    self.attrib = {}
    self.keys = self.attrib.keys()
    self.innerHTML = ''
    self.tag = ''
    self.start = -1
    self.end = -1

class Parser():
  def __init__(self,infile):
    self._root = open(infile).read()
    self._root = self._root[self._root.find('<body'):self._root.find('</body')]
    self._root = self._root[self._root.find('>')+1:]
    self.tags = set()
    i,j=0,0
    self.tagstarts = {}
    self.tagends = {}
    self.elements = []
    while i<len(self._root): 
      i = self._root.find('<',i)
      j = self._root.find('>',i)
      if not j>i: break
      tag = self._root[i+1:j].split()[0]     
      if tag.startswith('/'):
        tag = tag[1:]
        self.tagends[tag][-1-self.tagends[tag][::-1].index(None)]=((i,j))
      else:
        self.tags.add(tag)
        if tag in self.tagstarts:
          self.tagstarts[tag].append((i,j))
          self.tagends[tag].append(None)
        else:
          self.tagstarts[tag]=[(i,j)]
          self.tagends[tag]=[None]
      i=j+1
    self.getElements()
  def getElements(self):
    for tag in self.tags:
      if not tag in ['img']:
        for i in range(len(self.tagstarts[tag])):
          element = Element()
          element.startTag = self.tagstarts[tag][i][0]
          element.endTag = self.tagstarts[tag][i][1]
          tagData = self._root[element.startTag+1:element.endTag].replace("\"","")
          element.tag = tag
          element.attrib=dict([a.split('=') for a in tagData.split() if '=' in a])
          element.start = element.startTag
          element.end = self.tagends[tag][i][1]
          element.innerHTML = (self.tagstarts[tag][i][1]+1,self.tagends[tag][i][0])
          self.elements.append(element)
  def getText(self,tag,key,val):
    element = [element for element in self.elements if element.tag==tag and element.attrib[key]==val]
    if len(element):
      element = element[0]
      start,end = element.innerHTML
      text = self._root[start:end]
      i,j=0,0
      while i<len(text):
        i  = text.find('<',i)
        j = text.find('>',i)
        if i<0 or j<0: break
        text = text[:i]+text[j+1:]
      i,j=-1,-1   
      while i<len(text):
        i  = text.find('&',i+1)
        j = text.find(';',i+1)
        if i<0 or j<0: break
        if text[i+1:j] in ent.keys(): text = text[:i]+ent[text[i+1:j]]+text[j+1:]
    else:
      text=None 
    return text

Note the code does not check if the html code is correct. it also work only for the body part, and ignore img tags.
   

Howto export points to a Google Earth KML file

The problem:
A table of xyz data points can be viewed spatially with different methods:
GMT and ArcView are just a couple of options. I often get tables of spatial data points in Microsoft Excel (xp ver.) format or CSV of some kind. An easy (and free) way to visualize the data points on Windows would be to export the data and view it on Google Earth.
The solution:
points2kml (115kb).
In order to ease the work-flow I wrote a VBA macro add-in for Excel that exports selected fields (longitude, latitude, name and data) to a Google Earth KML file. the point2kml add-in and can be found on my scripts page.

Using hillshade image as intensity (improved matplotlib shade)

Matplotlib module enables hillshade method (shade) using a LightSource class (v 0.99).
The problem is it uses the data itself as intensity and data. It is very useful for viewing a DEM but sometimes you would like the DEM as intensity underlying some other data. Another problem is that the shade method is producing a very light colored image sometimes even white where intensity is high.
I used as an example a DEM derived from SRTM v4 data acquired at the International  Centre for Tropical  Agriculture (CIAT - http://srtm.csi.cgiar.org) the hillshade production was made using LightSource class with azimuth of 165 deg. and altitude of 45 deg.)

DEM - gist_earth color scheme

Hill-shade (azdeg-165,altdeg-45)

matplotlib shade method

My shade method

The difference in the shading colors derived from the method used to produce it. While the matplotlib method uses "hard light" method I use a "soft light" method. the matplotlib is converting the RGB colors to HSV and then calculate the new saturation and value according to the intensity. I use a formula based on the description of ImageMagick's pegtop_light.which is much faster as it is a single formula. Another advantage is the option to use a separate layer as the intensity and another as the data used for colors.
The modified functions are hillshade and set_shade as follows:

#!/bin/env python
from pylab import *
def set_shade(a,intensity=None,cmap=cm.jet,scale=10.0,azdeg=165.0,altdeg=45.0):
''' sets shading for data array based on intensity layer
  or the data's value itself.
inputs:
  a - a 2-d array or masked array
  intensity - a 2-d array of same size as a (no chack on that)
                    representing the intensity layer. if none is given
                    the data itself is used after getting the hillshade values
                    see hillshade for more details.
  cmap - a colormap (e.g matplotlib.colors.LinearSegmentedColormap
              instance)
  scale,azdeg,altdeg - parameters for hilshade function see there for
              more details
output:
  rgb - an rgb set of the Pegtop soft light composition of the data and 
           intensity can be used as input for imshow()
based on ImageMagick's Pegtop_light:
http://www.imagemagick.org/Usage/compose/#pegtoplight'''
  if intensity is None:
# hilshading the data
    intensity = hillshade(a,scale=10.0,azdeg=165.0,altdeg=45.0)
  else:
# or normalize the intensity
    intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
# get rgb of normalized data based on cmap
  rgb = cmap((a-a.min())/float(a.max()-a.min()))[:,:,:3]
# form an rgb eqvivalent of intensity
  d = intensity.repeat(3).reshape(rgb.shape)
# simulate illumination based on pegtop algorithm.
  rgb = 2*d*rgb+(rgb**2)*(1-2*d)
  return rgb

def hillshade(data,scale=10.0,azdeg=165.0,altdeg=45.0):
  ''' convert data to hillshade based on matplotlib.colors.LightSource class.
    input:
         data - a 2-d array of data
         scale - scaling value of the data. higher number = lower gradient
         azdeg - where the light comes from: 0 south ; 90 east ; 180 north ;
                      270 west
         altdeg - where the light comes from: 0 horison ; 90 zenith
    output: a 2-d array of normalized hilshade
'''
  # convert alt, az to radians
  az = azdeg*pi/180.0
  alt = altdeg*pi/180.0
  # gradient in x and y directions
  dx, dy = gradient(data/float(scale))
  slope = 0.5*pi - arctan(hypot(dx, dy))
  aspect = arctan2(dx, dy)
  intensity = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(-az - aspect - 0.5*pi)
  intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
  return intensity

Example of use:
One can save the code to a file named say: shading.py
now say we have a 4 byte float DEM data in a 560 lines 420 samples binary file. in a python code:

from pylab import *
from shading import set_shade
from shading import hillshade
dem = fromfile('DEM.dem',dtype=float32).reshape(560,420)
rgb = set_shade(dem,cmap=cm.gist_earth)
imshow(rgb)

will produce the "my shade method" image as above.

say we have a data to be plot using the DEM data as intensity:

replace the line before last with:
rgb = set_shade(data,intensity=hillshade(dem),cmap=cm.gist_earth)

howto edit multiple files via command line

In order to replace a regexp in multiple files using the command line use:

sed -i -e 's/regexp/new_regexp/' file_name

where file_name can be also a regexp (e.g. *.txt)

HowTo view a ROI PAC file? (or "replacement for mdx")

The Problem:
I can't get mdx/dgx in order to view my ROI PAC files.
The Solution:
a Python script called roiview (Get RoiView at ) wich is one of my other scripts.
This script uses pylab and can read, display and save images of ROI PAC quite similar to mdx/dgx.
it can handle *.unw, *.cor, *.slc, *.int and *.dem files.
just download, extract and use:
"roiview -h" for manual.
Update! Oct 13, 2010: A new version of RoiView (0.75) is now available at:

HowTo Create a SRTM DEM for ROI PAC Processing

When Processing interferograms, usually, you will also require a DEM. The NASA's SRTM project is freely available from different places in different formats and levels. The latest version for now (IV) is available from CGIAR Consortium for Spatial Information.
in order to ease the DEM production for a research, I wrote a small Python script. The script will create a DEM with 1 deg to each side of the center of scene requested. It will create a ROI PAC header as well. there is an option to create a JPG image of shaded relief of the imag. You will need Python to be installed and also PyLab and PIL modules.

simply download SRTM.py (zip,2.9kb) from my scripts page, change permission to execute, and run it...

Don't forget to credit CGIAR for the data. Citations should be made as follows: Jarvis A., H.I. Reuter, A. Nelson, E. Guevara, 2008, Hole-filled seamless SRTM data V4, International Centre for Tropical Agriculture (CIAT), available from http://srtm.csi.cgiar.org/

Howto view an image in Google Earth

The Problem:
I have an image I wish to view on Google Earth.

The Solution:
I assume you have an image in jpg or tiff or gif format. and you have a world file for it or can create one. Then, all you need to do is to run the image2kml.csh script.
just download (~1kb), apply execution privilege (chmod +x [ScriptFileName])
and run it. or visit my scripts page
Tip:
gif format image can support transparency you can use GIMP or other image manipulation program to turn unwanted parts of the image to be transparent so you can see the Google Earth layover under the image.

Whoto know if a long process is done?

The Problem:
I sometimes run a very long process which can even reach up to few days.
Instead of looking and checking if the process has done, I needed a solution to get a message from the process itself.
One way to do it is by using mail like this:
echo "some message" | mail -s "subject of mail" yourmail@wherever
Unfortunately, for some reason it doesn't work for me.
if you have any solution - let me know. I used a workaround...

The Solution (Option I):
Using Thunderbird and cshell and a tip from here (look at the comments), I built a script that can be added to a command line or another script which run's the process.
needed ingredients:
- Thunderbird
- xdotools

usage: simply copy the command below, edit where needed and run:

mail_script.csh [few message body words]

mail_script.csh:

#!/bin/csh
# this will send a mail with subject:"mission accomplished"
# and command line variables (argv) as the message body
# the script require Thunderbird and xdotool installed
# Just change the "to" field in the command below and
# set a program to run it at the end of process
# by Ran Novitsky Nof 2009 @ BGU
set body = `echo $argv`
thunderbird -compose "to='yourmail@wherever',subject='Mission Accomplished',body=$body" & ; sleep 5 ; set WID=`xdotool search --title
"Compose" head -1`; xdotool windowfocus $WID;xdotool key ctrl+Return

 The Solution (Option II):
  Option I needs Thunderbird to be installed and it uses gui so somtimes there are conflicts with user activity. an alternative solution is to install ssmtp (most Linux distros has it as a package). after installing edit /etc/ssmtp/ssmtp.conf file:

most importent lines to edit:
Mailhub - you smtp server (e.g smtp.gmail.com or any other smtp server)
AuthUser - see conf file link above
AuthPass - see conf file link above
AuthMethod - see conf file link above
FromLineOverride - remove the # at the start of line

edit /etc/ssmtp/revaliases:
add a line:
root:[user]@[mailaccount]:[smtpserver]:[port]

and run the script mail_script.csh:

# this will send a mail with subject:"mission accomplished"
# and command line varibles (argv) as the message body
# the script requier ssmtp installed
# Just change the -F , SendTo and "To" fields in the command below and set a program to run it at the end of process
#
# by Ran Novitsky Nof 2010 @ BGU
set body = `echo $argv`
sudo ssmtp -F"name of sender" recipient@mailserver << END
To:recipient@mailserver
Subject:Mission accomplished
$body
END

Tag Clouds

The Problem:
How to replace the tag list (aka labels) to a tag cloud?
The Solution:
Details in here you can see the outcome on the sidebar of my blog.
Tip:
To add a heading to the tags cloud simply add another <div> before the div of widget-content:

new div: <div><h2 class="'title'">Labels Cloud</h2></div> (in red is the title I selected for the gadget, can be anything actually)

original script (only the first few lines):<b:widget id="'Label1'" locked="'false'" title="'Labels'" type="'Label'"><b:includable id="'main'"> [new div in here] <div class="'widget-content'"><div id="'LabelDisplay'/"></div><script language="'javascript'" type="'text/javascript'">