Tuesday, March 8, 2016

Landsat QGIS tutorial -- Import and Band Ratio

This was written as a rough work flow for a colleague working on a project:

1) Create a personal account at http://earthexplorer.usgs.gov/ , mark your interest area and dates. Do not search to large date periods, only the first 100 results will be shown. Then choose Landsat-8 under Datasets>Landsat Archive:



2) Under Additional Criteria you can choose cloud cover threshold, but it seems to me that it is best to search for all scenes. A scene may have 80% cloud cover, but possibly just your glacier may be cloudfree on that image. There is very few cloud free images for Svalbard!

3)  In the result list, choose the scenes you want to download. Either add them to the bulk download, or use "download option" to choose a single scene. Get the GeoTiff file: 


4) The Landsat files are a Tif-file for each channel, but for displaying band combinations you need to merge these into one file:

5) Extract the Landsat zip-file and then find the merge dialog in QGIS: 


6) Choose all bands of the Landsat image to be merged. Be sure to have "layer stack" chosen!

7) For displaying a rgb image, you have to check which band is red, green or blue. This varies between Landsat sensors, consult this page for finding out. For Landsat 8, rgb is channel 4,3,2

Choose the property dialog (right-click the Landsat layer and choose "property" to view various band combination, in this example for Landsat-8 it is rgb, i.e. the regulart visual view by assigning 4,3,2 to r,g,b:


8) Band Ratio images: For discerning surfaced, often ratio between various bands are used. A known one is the NDVI index 
For glaciers, GLIMS suggests a few band ratio approaches, like "(TM3 / TM5)>2.6" Be aware that band numbers do not correspond to the same wavelength bands on different sensors, not even between different Landsat sensors! Check which actually wavelengths bands your formula operates with.

To create a band ratio image, use the QGIS raster calculator. This example calculates an output raster with (band4 / band5 )>2.6, but you could do any other operation you like

Reshaping Polygons in QGIS

We have shapefile giving the glacier outlines on Svalbard, now after a few years we want to update these outlines since glaciers have changed in extent:


For reshaping this polygon, activate the advanced digitizing menue in QGIS:


Then select the shapefile layer, choose "toggle editing" to make the shapefile layer editable and choose the "Reshape Feature" Button:

 Now start digitizing, first click outside the polygon, crossing the old glacier front, the digitize the new glacier front, and at the end again cross the old front line as indicated here. Below you see this for e retreating glacier, similarly an advancing glacier you cross the old front line before digitizing the new front. THe principle is nicely explained at this blog post. You can zoom into the image use the arrow keys to move the image while digitzing:

Start Digitizing:

End of Digitizing: 

For an advancing glacier you start and end by crossing the polygon like this, again look at this blog post to understand the principle: 



When you now right click, the polygon is adjusted to the new front line:

Now select the polygon you just edited, right click the layer and choose "open attribute table. In the table you choose to show the selected polygon, where you can update the relevant columns (year, date, analyst, source2 and sensor): 






 

Wednesday, February 24, 2016

Creating a Texture Image with a GLCM Co-Occurence Matrix using Scikit-Image and Python

Many clustering methods like the previously mentioned k-means cluster all the image values without considering any possible relation or patterns between neighbouring pixels. Texture is a measure considering this spatial relation --- an area with similar absolute values may visually be very different showing a different texture pattern (Source: Wikipedia)





The sci-kit image package for Python allows to calculate the Grey Level Co-Occurence Matrix (GCLM) and its texture measures as described on this page. However, this description, just as in the previous k-means example, leaves out the process of creating the texture image itself.

On this page you find a very good tutorial explaining GCLM and texture. A texture image is created by moving a window over you raster replacing the centre pixel covered by this window with the calculated texture measure of the moving window, as you see in this image taken from Myrka Hall-Beyers excellent tutorial:




A 7x7 window for any pixel would be defined by

glcm_window = sarraster[i-3: i+4, j-3 : j+4]


And the GCLM for a one-pixel offset [1] to the right [0] for this window as documented here would then be

glcm = greycomatrix(glcm_window, [1], [0],  symmetric = True, normed = True )


Now you have to calculate this for every single pixel. The only way for making this specific calculation I can come up with is the slow process of looping through each pixel. Since some software implementations (like the SNAP SAR Toolbox) make this much faster, there must be some other way, so maybe a reader of this post has an idea. (I asked this question also at http://stackoverflow.com/questions/35551249/implementing-glcm-texture-feature-with-scikit-image-and-python ).

The following is only practical for small raster images, named sarraster in this case:


for i in range(sarraster.shape[0] ):
    print i,
    for j in range(sarraster.shape[1] ):
        
        #windows needs to fit completely in image
        if i <3 or j <3:        
            continue
        if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4):
            continue
        
        #Define size of moving window
        glcm_window = sarraster[i-3: i+4, j-3 : j+4]
        #Calculate GLCM and textures
        glcm = greycomatrix(glcm_window, [1], [0],  symmetric = True, normed = True )

        # Calculate texture and write into raster where moving window is centered
        contrastraster[i,j]      = greycoprops(glcm, 'contrast')
        dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity')
        homogeneityraster[i,j]   = greycoprops(glcm, 'homogeneity')
        energyraster[i,j]        = greycoprops(glcm, 'energy')
        correlationraster[i,j]   = greycoprops(glcm, 'correlation')
        ASMraster[i,j]           = greycoprops(glcm, 'ASM')
        glcm = None
        glcm_window = None
        


The following is the original image with various texture images. This is just a first test example, so no clear patterns are seen in this image:



The complete script for creating the image above. For a 1100 x 1400 raster this script takes about an hour while in the SNAP SAR Toolbox such calculations take < 1 minute, so I wonder about alternative implementations:


# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import gdal, gdalconst
import numpy as np
from skimage.feature import greycomatrix, greycoprops

#Read SAR image into Numpy Array
filename = "//mnt//glaciology//SARimage.jpg"
sarfile = gdal.Open(filename, gdalconst.GA_ReadOnly)
sarraster = sarfile.ReadAsArray()


#Create rasters to receive texture and define filenames
contrastraster = np.copy(sarraster)
contrastraster[:] = 0

dissimilarityraster = np.copy(sarraster)
dissimilarityraster[:] = 0

homogeneityraster = np.copy(sarraster)
homogeneityraster[:] = 0

energyraster = np.copy(sarraster)
energyraster[:] = 0

correlationraster = np.copy(sarraster)
correlationraster[:] = 0

ASMraster = np.copy(sarraster)
ASMraster[:] = 0

# Create figure to receive results
fig = plt.figure()
fig.suptitle('GLCM Textures')

# In first subplot add original SAR image
ax = plt.subplot(241)
plt.axis('off')
ax.set_title('Original Image')
plt.imshow(sarraster, cmap = 'gray')


for i in range(sarraster.shape[0] ):
    print i,
    for j in range(sarraster.shape[1] ):
        
        # windows needs to fit completely in image
        if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4):
            continue
        
        # Define size of moving window
        glcm_window = sarraster[i-3: i+4, j-3 : j+4]
        # Calculate GLCM and textures
        glcm = greycomatrix(glcm_window, [1], [0],  symmetric = True, normed = True )

        # Calculate texture and write into raster where moving window is centered
        contrastraster[i,j]      = greycoprops(glcm, 'contrast')
        dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity')
        homogeneityraster[i,j]   = greycoprops(glcm, 'homogeneity')
        energyraster[i,j]        = greycoprops(glcm, 'energy')
        correlationraster[i,j]   = greycoprops(glcm, 'correlation')
        ASMraster[i,j]           = greycoprops(glcm, 'ASM')
        glcm = None
        glcm_window = None
        
        


texturelist = {1: 'contrast', 2: 'dissimilarity', 3: ' homogeneity', 4: 'energy', 5: 'correlation', 6: 'ASM'}
for key in texturelist:
    ax = plt.subplot(2,3,key)
    plt.axis('off')
    ax.set_title(texturelist[key])
    plt.imshow(eval(texturelist[key] + "raster"), cmap = 'gray')

plt.show()

Tuesday, February 23, 2016

K-Means Clustering of a Satellite Images using Scipy

Scipy can be used nicely to apply k-means classification on a satellite image. Most blog entries I find (such as this here) describe how to apply k-means to x/y data but not how to get a image classified into these k-means clusters.

In a first step you load the satellite image into a raster, but then flatten this raster to convert the row x column raster into a one-dimensional array

>>> flatsarraster = sarraster.flatten()
>>> sarraster.shape
(1486, 1109)
>>> flatsarraster.shape
(1647974,)
>>> sarraster
array([[-17.93579102, -17.26957893, -17.6971817 , ...,  -3.9816401 ,
         -4.78030491,  -5.49047613],
       [-15.8222208 , -18.36615753, -17.06410408, ...,  -5.16635323,
         -4.9837966 ,  -5.192698  ],
       [-13.59493637, -15.24811172, -13.07635021, ...,  -3.1499095 ,
         -3.1499095 ,  -5.5106535 ],
       ..., 
       [-15.59587097, -15.20659828, -14.20381069, ...,  -6.38638067,
         -6.44060993,  -3.72703171],
       [-15.73380089, -15.21032333, -14.21171951, ...,  -5.8094368 ,
         -6.5339241 ,  -4.01910543],
       [-15.65506649, -14.98183632, -14.98183632, ...,  -6.26982498,
         -5.8094368 ,  -4.90749454]], dtype=float32)
>>> flatsarraster
array([-17.93579102, -17.26957893, -17.6971817 , ...,  -6.26982498,
        -5.8094368 ,  -4.90749454], dtype=float32)

The command centroids, variance = kmeans(flatsarraster, [number of cluster]) then uses the flattened raster and calculates the centroids of the wanted numbers of clusters, in this case 8 clusters:
>>> from scipy.cluster.vq import *
>>> centroids, variance = kmeans(flatsarraster, 8)
>>> centroids
array([  0.44888815, -11.17466545, -20.26683044,  -5.95262289,
       -13.68685532,  -8.66477108,  -3.06361866, -16.28946304], dtype=float32)

 Now the command code, distance = vq(flatsarraster, centroids) gives a one dimensional array "code" of the same length as the flattened raster, so each position does not contain the number of the original image, but the number of the cluster it belongs to.

>>> flatsarraster
array([-17.93579102, -17.26957893, -17.6971817 , ...,  -6.26982498,
        -5.8094368 ,  -4.90749454], dtype=float32)
>>> code
array([7, 7, 7, ..., 3, 3, 3])
>>> 

So in a third step we reshape this array "code" into the "row x column" raster and get a image showing the clusters.

codeim = code.reshape(sarraster.shape[0], sarraster.shape[1])

Here the code in total, in this case creating a loop to classify 2 to 8 clusters:

# -*- coding: utf-8 -*-


import gdal, gdalconst
from scipy.cluster.vq import *
import matplotlib.pyplot as plt

# Read SAR image into Numpy Array
filename = "//mnt//glaciology//processedSARimages//SARimage.tif"
sarfile = gdal.Open(filename, gdalconst.GA_ReadOnly)
sarraster = sarfile.ReadAsArray()

# Flatten image to get line of values
flatsarraster = sarraster.flatten()

# Create figure to receive results
fig = plt.figure()
fig.suptitle('K-Means Classification')

# In first subplot add original SAR image
ax = plt.subplot(241)
plt.axis('off')
ax.set_title('Original Image')
plt.imshow(sarraster, cmap = 'gray')

# In remaining subplots add k-means classified images
for i in range(7):
    print "Calculate k-means with ", i+2, " cluster."
    
    #This scipy code classifies k-mean, code has same length as flattened
    #SAR raster and defines which class the SAR value corresponds to
    centroids, variance = kmeans(flatsarraster, i+2)
    code, distance = vq(flatsarraster, centroids)
    
    #Since code contains the classified values, reshape into SAR dimensions
    codeim = code.reshape(sarraster.shape[0], sarraster.shape[1])
    
    #Plot the subplot with (i+2)th k-means
    ax = plt.subplot(2,4,i+2)
    plt.axis('off')
    xlabel = str(i+2) , ' clusters'
    ax.set_title(xlabel)
    plt.imshow(codeim)
    
plt.show()

The resulting images, each color indicating a cluster:

Friday, November 20, 2015

Parsing a xml-file in Python with a Conditional Search to get Corner Coordinates of Satellite Images

I still find it somewhat tedious to search through xml-files. The following describes a solution on how to find variables in a xml-files that meet certain conditions:

A Radarsat satellite image comes with a product.xml file. In the depth of the xml-hierarchy, there is a long list of image tie points providing coordinates for selected image pixels:


From this long list, I want to retrieve those tie points being the corner coordinates of this image. In order to do this, I need to read from the xml-file

1) the number of pixels and the number of lines of this image in order to determine the corners
2) the tiepoints containing latitude and longitude information for these corners

A good starting point to understand how xml-files can be parsed in Python can be found here.

Using the elementtree module in Python, I first parse the xml file and get the root of the xml-file:

import xml.etree.ElementTree as etree      
tree = etree.parse(xml_file_extracted)
root = tree.getroot()

I can look at the immediate subelements of the root level using a simple loop:

>>> root = tree.getroot()
>>> for child in root:
...     print child.tag, child.attrib
... 
{http://www.rsi.ca/rs2/prod/xml/schemas}productId {}
{http://www.rsi.ca/rs2/prod/xml/schemas}documentIdentifier {}
{http://www.rsi.ca/rs2/prod/xml/schemas}sourceAttributes {}
{http://www.rsi.ca/rs2/prod/xml/schemas}imageGenerationParameters {}
{http://www.rsi.ca/rs2/prod/xml/schemas}imageAttributes {}
>>> 

Note that the name of the element is preceded by "{http://www.rsi.ca/rs2/prod/xml/schemas}", so if you search an element named rasterAttribute, the full name is in fact "{http://www.rsi.ca/rs2/prod/xml/schemas}rasterAttribute"

I now need to search for this element "rasterAttribute." I do not know its exact path within the xml-file, but I know from looking at the file that this element contains the number of pixels and the number of lines. Using "root.iter" I can iterate through all elements with a given name, I use it here as a simple way to find the one rasterAttribute element:


    for attributes in root.iter('{http://www.rsi.ca/rs2/prod/xml/schemas}rasterAttributes'):
        attribute_list = attributes.getchildren()


The attribute_list gives me all the subelements, the children of the rasterAttribute element:

[<Element '{http://www.rsi.ca/rs2/prod/xml/schemas}dataType' at 0x6158588>,
<Element '{http://www.rsi.ca/rs2/prod/xml/schemas}bitsPerSample' at 0x61585f8>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}numberOfSamplesPerLine' at 0x6158630>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}numberOfLines' at 0x6158668>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}sampledPixelSpacing' at 0x61586a0>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}sampledLineSpacing' at 0x61586d8>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}lineTimeOrdering' at 0x6158710>,
 <Element '{http://www.rsi.ca/rs2/prod/xml/schemas}pixelTimeOrdering' at 0x6158748>]

If I look at the rasterAttribute element in the xml-file, it would look like this:



It is the  numberOfSamplesPerLine and the numberOfLines elements I want to retrieve. I could simply get these from the attribute_list using:

attribute_list[2].text
attribute_list[3].text

attribute_list[2] reads the third element of the list being "numberOfSamplesPerLine" and the ".text" is the ElementTrees command to access its valus.

This approach of accessing the attribute_list would work if the structure of the xml-file always remained the same. However, different products from Radarsat have slightly different xml-files. For SLC products, the wanted elements are for example at the fourth and fifth position of the attribute_list.  The two lines above will therefore not always retrieve the correct element, so one needs to check that the correct element is read. In order to do this, I loop through the attribute_list, which contains all subelements within the rasterAttribute element, and if the element has the correct name, then I read the value into a variable:


for attribute in attribute_list:
      if attribute.tag == "{http://www.rsi.ca/rs2/prod/xml/schemas}numberOfSamplesPerLine":
           numberOfSamplesPerLine = float(attribute.text) - 1
      if attribute.tag == "{http://www.rsi.ca/rs2/prod/xml/schemas}numberOfLines":
           numberOfLines = float(attribute.text) - 1

Why the minus 1? The first pixel at line 1 and pixel 1 has the location (0,0), so the number of lines and rows is one higher than the "address" of the pixel.

Now I know the number of pixels and the number of lines of my image. The next step is to search through all elements called "imageTiePoint", finding the ones matching the corners of the image. The upper left corner is then the one where the pixel location is at line 0 and at pixel 0. The upper right at line 0 and pixel "numberOfSamplesPerLine" and so on ... Again, I loop through all elements called imageTiePoint, create a list "child_list" containing all subelements and then read the pixel and line values. If these match with the corners, I read the corresponding latitude and longitude:


for attribute in attribute_list:
    for tiepoint in root.iter('{http://www.rsi.ca/rs2/prod/xml/schemas}imageTiePoint'):
        child_list = tiepoint.getchildren()
        
        line      = float(child_list[0][0].text)
        pixel     = float(child_list[0][1].text)
        latitude  = float(child_list[1][0].text)
        longitude = float(child_list[1][1].text)
        
        
        #check if upper left
        if (line == 0.0) and (pixel == 0.0):
            sar_upperleft_x = longitude
            sar_upperleft_y = latitude
        
        #check if upper right
        if (line == 0.0) and (pixel == numberOfSamplesPerLine):
            sar_upperright_x = longitude
            sar_upperright_y = latitude
            
                
        #check if lower left    
        if (line == numberOfLines) and (pixel == 0.0):
            sar_lowerleft_x = longitude
            sar_lowerleft_y = latitude
            
                
        #check if lower right
        if (line == numberOfLines) and (pixel == numberOfSamplesPerLine):
            sar_lowerright_x = longitude
            sar_lowerright_y = latitude


For these I actually do access the values through its position within the child_list since within the imageTiePoint element they appear to be the same for all Radarsat products. It would be safer of course to check the name of the element the same way as I did for the rasterAttribute element.

I can with these values create a WKT representation of a polygon, which is then the bounding box of the satellite image:
    wkt_sarimage = "POLYGON((" + str(sar_upperleft_x) + " " +  str(sar_upperleft_y) + "," \
                  + str(sar_upperleft_x) + " " + str(sar_lowerright_y) + \
                  "," + str(sar_lowerright_x) + " " + str(sar_lowerright_y) + "," + str(sar_lowerright_x) \
                  + " " + str(sar_upperleft_y) + "," + str(sar_upperleft_x) + " " +  str(sar_upperleft_y) + "))"

Friday, October 23, 2015

Batch Downloading Sentinel Images from ESA's Scihub with Python

ESA's Sentinel Satellite Images can be downloaded from https://scihub.esa.int/ with batch scripting as described here: https://scihub.esa.int/userguide/5APIsAndBatchScripting

The following describes a simple script on how to batch-download Sentinel images over an area of interest using Python. The search request in this version below is hard-coded into the script, but solves the main-task of downloading these files. The following is a first attempt and certainly will be improved. I will update this page then.

Please note that scihub has now two alternative download locations. The one called "dhus" contains for now data a longer time back, but is not very stable. The one called "apihub" contains the last 30 days, but is dedicated for script download and more stable. For apihub replace "dhus" in all links below.


Authenticate at Webpage 

Since scihub needs registration, any web-request passed to scihub needs authentication. So first you have to open your access to this page. The tools within Python's urllib library needed for this are described at this page.

    
    # The lower level url to authenticate
    'https://scihub.copernicus.eu/dhus/'
    
    # Create a password manager
    # as described at 
    # https://docs.python.org/2.7/howto/urllib2.html#openers-and-handlers
    username = 'username'
    password = 'password'
    p = urllib2.HTTPPasswordMgrWithDefaultRealm() 
    p.add_password(None, url, username, password)
    handler = urllib2.HTTPBasicAuthHandler(p)

    # create "opener" (OpenerDirector instance)
    opener = urllib2.build_opener(handler)
    
    # Install the opener.
    # Now all calls to urllib2.urlopen use our opener.
    urllib2.install_opener(opener)


Create the query and save resulting xml-response

They query itself is a weblink created as described here: https://scihub.esa.int/twiki/do/view/SciHubUserGuide/5APIsAndBatchScripting#Open_Search

As an example, the following weblink searches for all GRD type products ingested in scihub between the given dates within the area of interest defined by the polygon:


The first issue is, however, to reformat the link above in order for te urllib library to be able to read it. Pasting the link above into the browser will give you the list of results. However, inspecting the link closely shows that what is actually passed on is a link, where some but not all special characters are replaced: 

In order to replace the needed characters, you can use the urllib.quote(string[, safe]) command. The safe parameter defines all the characters, which are not to be replaced, and after some trial and error I find that this string need to be
':()[]/?=,&'

For the url therefore you use urllib.quote, then read the result and write it to an xml file:


    urlrequest = urllib.quote('https://scihub.esa.int/dhus/search?q=productType:GRD AND ingestionDate:[2015-09-30T00:00:00.000Z TO 2015-09-30T23:59:59.999Z] AND footprint:"Intersects(POLYGON((-16.11154750824 69.190260085686,-2.7521725082398 69.190260085686,-2.7521725082398 76.370271635301,-16.11154750824 76.370271635301,-16.11154750824 69.190260085686)))"&rows=10000&start=0', ':()[]/?=,&')
    
    # Read page response into variable
    page = urllib2.urlopen(urlrequest).read()
    
    # Write returned xml list to textfile
    textfile = open('/home/max/Documents/test.xml', 'w')
    textfile.write(page)
    textfile.close()

This result creates an xml-file containing a list of matching products.

Parsing the xml file

The next step is now to parse this xml file for the resulting and matching files and download these. How to do this, I found nicely described in "Dive Into Python" It is worth playing a bit with the different examples from chapter 12 to understand how the ElementTree library works.

Basically, you create a list containing all elements called "entry" since these contain the matching scenes. For each entry tag you need to find the Identifier in the "id" element and the title.

The id element may, for example contain the uuid "18f7993d-eae1-4f7f-9d81-d7cf19c18378", from which we then know that the file can be downloaded under the adress
https://scihub.copernicus.eu/dhus/odata/v1/Products('18f7993d-eae1-4f7f-9d81-d7cf19c18378')/$value

The name of the image itself is in the "title" element and is needed to name the downloaded Sentinel-file.

This code loops through all entries and finds the uuid element, the name of the Sentinel image and creates the link where the file can be downloaded from:

    tree = etree.parse('/home/max/Documents/test.xml')
    entries = tree.findall('{http://www.w3.org/2005/Atom}entry')
    for entry in range(len(entries)):
        uuid_element = entries[entry].find('{http://www.w3.org/2005/Atom}id')
        sentinel_link = "https://scihub.copernicus.eu/dhus/odata/v1/Products('" + uuid_element.text + "')/$value"
        title_element = entries[entry].find('{http://www.w3.org/2005/Atom}title')

Downloading each matching image


The last task is now to download each file. I found some hints on this here and adding this to the for-loop above, I get this code:

    tree = etree.parse('/home/max/Documents/test.xml')
    entries = tree.findall('{http://www.w3.org/2005/Atom}entry')
    for entry in range(len(entries)):
        uuid_element = entries[entry].find('{http://www.w3.org/2005/Atom}id')
        sentinel_link = "https://scihub.copernicus.eu/dhus/odata/v1/Products('" + uuid_element.text + "')/$value"
        title_element = entries[entry].find('{http://www.w3.org/2005/Atom}title')

        # Full path where downloaded file is to be stored
        destinationpath =  "/home/max/Documents/SentinelDownloads/" +    title_element.text + '.zip'
        
        #Open the file and read the data
        print "Downloading ", title_element.text
        downloadfile = urllib2.urlopen(sentinel_link)
        data =  downloadfile.read()
        #Write the data to the destination zipfile
        with open(destinationpath, "wb") as code:
           code.write(data)

Using all this now, you have to create your own query and replace the urlrequest variable. Add your username and password, as well as define your destinationpath:

The complete script


import urllib2, urllib
import xml.etree.ElementTree as etree

# Hard coded link returning list of matching scenes
# https://scihub.esa.int/twiki/do/view/SciHubUserGuide/5APIsAndBatchScripting#Open_Search

# Authenticate at scihub webpage
url =  'https://scihub.copernicus.eu/dhus/'
username = 'your_username'
password = 'your_password'
password_mgr = urllib2.HTTPPasswordMgrWithDefaultRealm()

password_mgr.add_password(None, url, username, password)
handler = urllib2.HTTPBasicAuthHandler(password_mgr)
opener = urllib2.build_opener(handler)
urllib2.install_opener(opener)

# The request query
urlrequest = urllib.quote('https://scihub.copernicus.eu/dhus/search?q=productType:GRD AND ingestionDate:[2015-09-30T00:00:00.000Z TO 2015-09-30T23:59:59.999Z] AND footprint:"Intersects(POLYGON((-16.11154750824 69.190260085686,-2.7521725082398 69.190260085686,-2.7521725082398 76.370271635301,-16.11154750824 76.370271635301,-16.11154750824 69.190260085686)))"&rows=10000&start=0', ':()[]/?=,&')
# Read the response into page and write it to a xml-file
page = urllib2.urlopen(urlrequest).read()
textfile = open('/home/max/Documents/test.xml', 'w')
textfile.write(page)
textfile.close()

#Parsing the xml-file, the entry tag contains the results
tree = etree.parse('/home/max/Documents/test.xml')
entries = tree.findall('{http://www.w3.org/2005/Atom}entry')
print "Number of Scenes Found: ", len(entries)
for entry in range(len(entries)):
    #The uuid element allows to create the path to the file
    uuid_element = entries[entry].find('{http://www.w3.org/2005/Atom}id')
    sentinel_link = "https://scihub.copernicus.eu/dhus/odata/v1/Products('" + uuid_element.text + "')/$value"
    
    #the title element contains the corresponding file name
    title_element = entries[entry].find('{http://www.w3.org/2005/Atom}title')
    
    #Destinationpath with filename where download to be stored
    destinationpath =  "/home/max/Documents/SentinelDownloads/" +    title_element.text + '.zip'
    
    print
    print "Scene ", entry + 1 , "of ", len(entries)
    
    #Check if file already was downloaded
    print sentinel_link
    if os.path.exists(destinationpath):
        print title_element.text, ' already downloaded'
        
        continue
    
    #Download file and read
    print "Downloading ", title_element.text
    downloadfile = urllib2.urlopen(sentinel_link)
    data =  downloadfile.read()
    
    # Write the download into the destination zipfile
    with open(destinationpath, "wb") as code:
        code.write(data)

print "Done downloading"

Thursday, August 28, 2014

Fixing a Map Projection Issue for the SNAP S1TBX Toolbox

[This post was written for the Nest SAR Toolbox, which is now replaced by the SNAP SAR Toolbox]

In the SNAP Sar Toolbox you have a choice of coordinate systems:

For some of these, however, SNAP reports the error "Axis Too Complex".

I compared a map projection that works with one where “Axis too complex” error is reported. Map projections that work contain these lines:
 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH]

While the non-working ones contain these lines, which apparently are wrong or unprocessable by SNAP
  AXIS["Easting", "North along 90 deg East"], 
  AXIS["Northing", "North along 0 deg"],



If you use the graph builder or command line (gpt) where the parametres are defined in an XML-file, you can manually replace these lines in the graph file (xml file) for the projection definition, then it works! An example for such an XML-file used by SNAP is here.

Another example, this one works

 PROJCS["WGS 84 / Australian Antarctic Lambert", 
  GEOGCS["WGS 84", 
    DATUM["World Geodetic System 1984", 
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
      AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH], 
    AUTHORITY["EPSG","4326"]], 
  PROJECTION["Lambert_Conformal_Conic_2SP", AUTHORITY["EPSG","9802"]], 
  PARAMETER["central_meridian", 70.0], 
  PARAMETER["latitude_of_origin", -50.0], 
  PARAMETER["standard_parallel_1", -68.5], 
  PARAMETER["false_easting", 6000000.0], 
  PARAMETER["false_northing", 6000000.0], 
  PARAMETER["scale_factor", 1.0], 
  PARAMETER["standard_parallel_2", -74.5], 
  UNIT["m", 1.0], 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH], 
  AUTHORITY["EPSG","3033"]]


EPSG3031 works not: 
  PROJCS["WGS 84 / Antarctic Polar Stereographic", 
  GEOGCS["WGS 84", 
    DATUM["World Geodetic System 1984", 
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
      AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH], 
    AUTHORITY["EPSG","4326"]], 
  PROJECTION["Polar Stereographic (variant B)", AUTHORITY["EPSG","9829"]], 
  PARAMETER["central_meridian", 0.0], 
  PARAMETER["Standard_Parallel_1", -71.0], 
  PARAMETER["false_easting", 0.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["Easting", "North along 90 deg East"], 
  AXIS["Northing", "North along 0 deg"], 
  AUTHORITY["EPSG","3031"]]