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"]]

Wednesday, July 23, 2014

Reading AMSR-2 Data into GeoTIFF raster using gdal_grid

The passive microwave data from the Advanced Microwave Scanning Radiometer (AMSR), available at https://gcom-w1.jaxa.jp/auth.html, is delivered in orbital swaths. The raster band contain the image data in various frequency range, and the geolocation information for 89GHz is stored in one containing latitude values for each pixel, and another channel containing longitude values for each pixel. This data needs processing before having a regular geocoded raster.

The following is not an exhaustive description, but extended notes on how to read AMSR-2 files, which possibly may be of   help to others trying to solve a similar task.

For some, the final commented script at our github page may be enough help, the text below tries to explain the most important steps in this script:

The various frequency bands of AMSR-2 have different resolution (see the product specs and the user manual ), we choose the L1R dataset, where the data is already processed to match the geolocation stored in the lat/long band for 89GHz.

Information on the hdf image can be retrieved with gdalinfo:

gdalinfo C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5

The metadata printed out contains the names of the subchannels, as seen here:
and information about a specific channel can be retrieved by using gdalinfo on exactly this subchannel-name from the metadata:

gdalinfo HDF5:"C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5"://Brightness_Temperature_(89.0GHz-A,V)

Opening one of the bands in QGIS, the data looks like this (all images in this post are © JAXA EORC ):


In this example we want to convert the data into the NSIDC sea ice raster (EPSG:3411) covering the Arctic areas, using the gdal_grid utility, which creates a regular grid out of a point collection.

Important to know: The 89GHz channel is divided into 89A and 89B, and only both together give the full resolution of about 5km. Each 89 GHz channel has an associated latitude and a longitude raster giving the coordinates for each pixel. For the L1R product, the geolocation information of other frequencies with lower resolution can be derived from the 89A longitude and latitude raster by using only their odd columns.

In the first step, I unzip the gz - zipfiles, then open the hdf file (*.h5). The various frequency bands are stored in Subdatasets, so you open the hdf file with gdal.open(), but then use again gdal.open() for a subdataset (for information on the bands you can run gdalinfo on the *.h5 files:


HDFfile = gdal.Open( r'C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5' )

HDF_bands = HDFfile.GetSubDatasets()

#HDF Subdatasets are opened just as files are opened:
HDF_SubDataset = gdal.Open(HDF_bands[channel][0])
HDF_SubDataset_array = HDF_SubDataset.ReadAsArray()

HDF_Lat89A = gdal.Open(HDF_bands[46][0])
HDF_Lat89A_array = HDF_Lat89A.ReadAsArray()
 
HDF_Lon89A = gdal.Open(HDF_bands[48][0])
HDF_Lon89A_array = HDF_Lon89A.ReadAsArray()


In the next step, a create a comma-separated file containing longitude, latitude, brightness values for each raster point. This comma separated file is then the input for gdal_grid. I loop through each pixel and write the three values (  longitude, latitude, brightness ) into a csv-file.

So for the 89GHz channel, I write both 89A and 89B to a csv-file:
 
#Add header line to textfile
textfile = open( AMSRcsv, 'w')
textfile.write('lon,lat,brightness\n')

## Loop through each pixel and write lon/lat/brightness to csv file
for i in range(rows):
    for j in range(cols):
        wgs84=pyproj.Proj("+init=EPSG:4326")
        EPSG3411=pyproj.Proj("+init=EPSG:3411")
     
        lonA = HDF_Lon89A_array[i,j]
        latA = HDF_Lat89A_array[i,j]

        # lon/lat written to file already projected to EPSG:3411
        (lonA_3411, latA_3411) = pyproj.transform(wgs84, EPSG3411, lonA, latA)
        brightnessA = HDF_Br89AH_array[i,j]* 0.01 #APPLYING SCALING FACTOR!

        lonB = HDF_Lon89B_array[i,j]
        latB = HDF_Lat89B_array[i,j]

        # lon/lat written to file already projected to EPSG:3411
        (lonB_3411, latB_3411) = pyproj.transform(wgs84, EPSG3411, lonB, latB)
        brightnessB = HDF_Br89BH_array[i,j]* 0.01 #APPLYING SCALING FACTOR!

        if 35 < latA < 90:
            textfile.write(str(lonA_3411) + ',' + str(latA_3411) + ',' + str(brightnessA) + '\n')

        if 35 < latB < 90:
            textfile.write(str(lonB_3411) + ',' + str(latB_3411) + ',' + str(brightnessB) + '\n')
       
textfile.close()


For the lower resolution channels, I use the odd numbers of the 89A long and lat channel:


 
#Add header line to textfile
textfile = open( AMSRcsv, 'w')
textfile.write('lon,lat,brightness\n')

## Loop through each pixel and write lon/lat/brightness to csv file
for i in range(rows):
    for j in range(cols):
        wgs84=pyproj.Proj("+init=EPSG:4326")
        EPSG3411=pyproj.Proj("+init=EPSG:3411")

        #For low resolution the odd columns of Lon/Lat89 array to be taken!
        lonA = HDF_Lon89A_array[(i) ,(j*2+1)]
        latA = HDF_Lat89A_array[(i) ,(j*2+1)]

        # lon/lat written to file already projected to EPSG:3411
        (lonA_3411, latA_3411) = pyproj.transform(wgs84, EPSG3411, lonA, latA)
        brightnessA = HDF_SubDataset_array[i,j]* 0.01 #APPLYING SCALING FACTOR!

        if 35 < latA < 90:
            textfile.write(str(lonA_3411) + ',' + str(latA_3411) + ',' + str(brightnessA) + '\n')

textfile.close()

Now I can almost run the gdal_grid, but as described at http://www.gdal.org/gdal_grid.html I need to create a xml file describing my comma-separated csv file.
 
<OGRVRTDataSource>
 <OGRVRTLayer name="GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H">
  <SrcDataSource>G:\AMSR\GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H.csv</SrcDataSource>
  <GeometryType>wkbPoint</GeometryType>
  <GeometryField encoding="PointFromColumns" x="lon" y="lat" z="brightness" />
 </OGRVRTLayer>
</OGRVRTDataSource>

This xml file above can be created in a python script in the following manner, some more info here:
root = ET.Element("OGRVRTDataSource")

OGRVRTLayer  = ET.SubElement(root, "OGRVRTLayer")
OGRVRTLayer.set("name", AMSRcsv_shortname)

SrcDataSource = ET.SubElement(OGRVRTLayer, "SrcDataSource")
SrcDataSource.text = AMSRcsv

GeometryType = ET.SubElement(OGRVRTLayer, "GeometryType")
GeometryType.text = "wkbPoint"

GeometryField = ET.SubElement(OGRVRTLayer,"GeometryField")
GeometryField.set("encoding", "PointFromColumns")

GeometryField.set("x", "lon")
GeometryField.set("y", "lat")
GeometryField.set("z", "brightness")

tree = ET.ElementTree(root)
tree.write(AMSRcsv_vrt)


Now we finally can run gdal_grid, either command line:

gdal_grid -a_srs EPSG:3411 -a average:radius1=4000:radius2=4000:min_points=1 -txe -3850000 3750000 -tye -5350000 5850000 -outsize 760 1120 -l GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H.vrt GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H.tif

or called from a Python script:

AMSRcsv_shortname =  GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H
AMSRcvs_vrt = GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H.vrt
AMSR_tif = GW1AM2_201301010834_032D_L1SGRTBR_1110110_channel89H.tif  

radius1 = 4000  
radius2 = 4000  
os.system('gdal_grid -a_srs EPSG:3411 -a average:radius1=' + str(radius1) + ':radius2=' + str(radius2) + ':min_points=1 -txe -3850000 3750000 -tye -5350000 5850000 -outsize 760 1120 -l ' + AMSRcsv_shortname + ' '  + AMSRcsv_vrt + ' ' + AMSR_tif)

The radius indicates in which distance around a given output raster point the algorithm searches for points falling into the NSIDC raster -- if too small, it will result in empty pixels, if too large there will be too much smoothing since many pixels are averaged into one.

The result is, finally, the part of the swath falling into the NSIDC-raster:


In a final step, I take all of such swaths for one day and average them into a full image of that given day, see the Average Daily function in the script for details (images  © JAXA EORC ).


One issue using gdal_grid is its very low performance regarding speed (see this comment ), one 89GHz band takes 10 minutes and a lower resolution band 2 minutes calculation time. This is then about 25 minutes for all channels of one hdf file, and since every day has about 20 files, this means 8 hours for one averaged daily raster. gdal_grid may therefore not always be feasible until the speed issue is improved.

Thursday, June 26, 2014

Reading and Map projecting raster data with geolocation arrays using gdal

The following describes how to geocode a raster band, whose the geolocation information is stored in one channel containing latitude values for each pixel, and another channel containing longitude values for each pixel. Opening the image, you see the unprocessed orbital swath.


Information on the hdf image can be retrieved with gdalinfo:

gdalinfo C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5

The metadata printed out contains the names of the subchannels, as seen here:
and information about a specific channel can be retrieved by using gdalinfo on exactly this subchannel-name from the metadata:

gdalinfo HDF5:"C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5"://Brightness_Temperature_(89.0GHz-A,V)

The tricky part is now, to transfer such raster bands into a map projection.

gdalwarp can  read separate latitude and longitude rasters, which is called geolocation arrays in gdal. the comment on this page brought me a step further although I think that the vrt-code and the commands given are not completely correct. One has to create a vrt file describing the image raster, also referring to the latitude and longitude bands. Then in addition a vrt file describing the latitude band and a vrt file describing the longitude band.

So I create a vrt-file (see here for info) named GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt:


   
     C:\Users\max\Documents\lon.vrt
     1
     C:\Users\max\Documents\lat.vrt
     1
     0
     0
     1
     1
    
   
   
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Brightness_Temperature_(89.0GHz-A,H)
      1
      
      
      
    
  




Then a file named lon.vrt :


  
    
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Longitude_of_Observation_Point_for_89A
      1
      
      
      
    
  



Finally a file named lat.vrt :

  
    
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Latitude_of_Observation_Point_for_89A
      1
      
      
      
    
  



With these files created and placed in the same directory, I can use gdalwarp, but note that you are calling the vrt file rather than the hdf file!
gdalwarp -geoloc -t_srs EPSG:4326 C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt C:\Users\max\Documents\test100.tif


The result is the following swath:


In order to get it onto the polar stereographic projection, in my case the NSIDC EPSG:3411 projection, I first subset the image:


gdal_translate -projwin -94 90 40 35 C:\Users\max\Documents\test100.tif C:\Users\max\Documents\test101.tif


which results in this image:



now I reproject into EPSG:3411:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3411 C:\Users\max\Documents\test101.tif C:\Users\max\Documents\test102.tif

resulting in this image:

Monday, June 16, 2014

Converting Coordinates between map projections using Python

Python can be used very nicely to quickly convert coordinates from one projection into another using the pyproj package, as usual available here. This example does it from command line, but you can built it into your script, of course, as well. The easiest way is to find the EPSG code for your projections, for example on this page. You then initialize the projections you want, in this example unprojected lat/long into EPSG:3575:

>>> EPSG4326 = pyproj.Proj("+init=EPSG:4326")
>>> EPSG3575 = pyproj.Proj("+init=EPSG:3575")


Alternatively you can use the Proj.4 definition, so the following two lines are doing the same, you can use either of them
EPSG3995 = pyproj.Proj("+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ")
EPSG3995 = pyproj.Proj("+init=EPSG:3995")


Now you can define the input coordinates and transform them, in this example you transform values from EPSG:4326 to EPSG:3575:

>>> lat = 79
>>> long = -5
>>> x, y = pyproj.transform(EPSG4326, EPSG3575, long, lat)
>>> x,y
(-317466.6103301885, -1184801.519458934)

Official documentation can be found here

Thursday, March 20, 2014

Quantarctica: An Antarctic GIS package

[This article was first contributed to the Case Studies section at the QGIS website in August 2013. Putting it here as well, as another advertisement!]

Quantarctica is a collection of Antarctic geographical datasets, such as base maps, satellite imagery, glaciology and geophysics data from data centres around the world, prepared for viewing in Quantum GIS (QGIS). The package is developed by the Norwegian Polar Institute, as a tool for the research community, for classrooms and for operational use in Antarctica – freely available for non-commercial purposes.


About the project



Screenshot from Quantarctica, showing one of the subglacial lakes datasets.

Quantarctica (Quantum GIS + Antarctica) was first developed for in-house use at the Norwegian Polar Institute, as a tool for our glaciologists. There was a need for a low or no cost complete GIS with essential datasets – ready-to-use, easy-to-use, functionality rich and with offline capabilities. QGIS seemed to be a perfect choice of GIS for the collected datasets.

Quantarctica has been used to examine geographical data from continental to local scales, for viewing scientific project data on top of base maps or with other scientific datasets, and to prepare maps for publications and proposals. Quantarctica has so far proven to be a great tool, and a very good alternative and supplement to other software used by the researchers. It has provided new opportunities for our researchers in their daily work.


Quantarctica is also useful when navigating on the Antarctic ice shelves thanks to the GPS tracking capabilities within QGIS. (Image by Peter Leopold)

Since Quantarctica first came in use by our glaciologists three years ago, there has been many requests in the research community outside the institute to share this product, and we started to develop a public and improved version to replace the in-house version. Following Antarctic field testing, and adding new relevant datasets, Quantarctica version 1.0 was finally completed and made available for download in July 2013.

Quantarctica is to be all about community effort. With contributions we aim to expand with data from other disciplines, such as oceanography, atmospheric sciences, geology and biology, and hope and believe that this tool can be useful for the Antarctic community – as a complete Antarctic GIS package.


Links


Quantarctica website: http://quantarctica.org/
Norwegian Polar Institute: http://www.npolar.no/en/


Monday, March 3, 2014

Geocoding NSIDC Sea Ice Concentration Maps

NSIDC sea ice projection maps have their own EPSG code for their map projection, which for the Artic is EPSG:3411. This map projection is part of the gdal libraries and also available in Quantum GIS.

The files themselves are delivered as flat binary strings and are converted to GeoTIFF files as described in a previous post . It took me a while to figure out how to apply the map projection to these files, so here is a short description on what to do:

Information on the NSIDC specific polar stereographic projection can be found on this web page. A figure from this page shows the extent:

Applying the map projection to my raster is not that difficult:


    driver = gdal.GetDriverByName("GTiff")
    outraster = driver.Create(outfile,  width, height,1, gdal.GDT_Int16 )
    
    spatialRef = osr.SpatialReference()
    spatialRef.ImportFromProj4('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs')
    outraster.SetProjection(spatialRef.ExportToWkt() )

An important note -- for some unknown reasons, the code does not work using spatialRef.ImportFromEPSG(3411) but it does with using the proj4 import as above, see also discussion here.

But somehow it took me a while to figure out how to set the geotransform (which is in a way the gdal version of the ESRI world file). But the NSIDC web page has all the necessary information in Table 2, where the corner coordinates are defined. The geotransform is defined as:



geotransform = (upperleft_x, pixelsize, rotation, upperleft_y, pixelsize, rotation)


which in our case -- using the values from Table 2, converted to metres and using negative directon for resolution --- becomes:
    driver = gdal.GetDriverByName("GTiff")
    geotransform = (-3850000, 25000 ,0 ,5850000, 0, -25000)
    outraster.SetGeoTransform(geotransform)

The resulting image shows the correct coordinates in lat/long when using the QuantumGIS "Coordinate Capture" Tool. The first picture shows the map projected file, for the picture below I use gdalwarp to project the map to EPSG:3575 which I prefer to use for the Barents Sea and Fram Strait:




For the reprojection from EPSG:3411 to EPSG:3575 I use this small function in my script -- it forces the output to keep the 25km resolution:
 
def EPSG3411_2_EPSG3575(infile):
    '''
    reprojects the infile from NSIDC 3411 to EPSG 3575
    outputfiel has 25km resolution
    '''
    
    (infilepath, infilename) = os.path.split(infile)
    (infileshortname, extension) = os.path.splitext(infilename)
    outfile = infilepath + '\\EPSG3575\\' + infileshortname + '_EPSG3575.tif'
    print ' Reproject ', infile, ' to ', outfile 
    os.system('gdalwarp -s_srs EPSG:3411 -tr 25000 -25000 -t_srs EPSG:3575 -of GTiff ' + infile + ' ' + outfile)


    
 

A close-up of Svalbard, also showing maximum and minimum of sea ice in a given month:

Thursday, February 27, 2014

Reading binary data (NSIDC Sea Ice Concentrations) into GeoTIFF raster with Python

I need to read NSIDC sea ice concentration charts from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data into a raster file. These files are according to the documentation "scaled, unsigned flat binary with one byte per pixel, and therefore have no byte order, or endianness. Data are stored as one-byte integers representing scaled sea ice concentration values.[...] The file format consists of a 300-byte descriptive header followed by a two-dimensional array of one-byte values containing the data."

The blog entry, which helped med finding out how to do this is this one but it had to be modified for my files.


Part One uses the Python struct command to unpack he binary file into a tuple
import struct, numpy, gdal

#Dimensions from https://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
height = 448
width = 304

#for this code, inspiration found at https://stevendkay.wordpress.com/category/python/
icefile = open(r"U:\SSMI\nt_20120101_f17_v01_n.bin", "rb")
contents = icefile.read()
icefile.close()

# unpack binary data into a flat tuple z
s="%dB" % (int(width*height),)
z=struct.unpack_from(s, contents, offset = 300)

nsidc = numpy.array(z).reshape((448,304))
nsidc = numpy.rot90(nsidc, 1)
The "unpack_from" considers the offset of the header bytes and the tuple "z" contains now all the values in one long row of values -- this one-dimensional string of values is converted into an array and the reshaped according to the dimensions of the raster. Now I have my values in a two-dimensional array, the next part of the code writes it into a GeoTIFF file:
#write the data to a Geotiff
driver = gdal.GetDriverByName("GTiff")
outraster = driver.Create('C:\Users\max\Desktop\\test4.tif', height, width, 1, gdal.GDT_Int16 )
raster = numpy.zeros((width, height), numpy.float) 
outraster.GetRasterBand(1).WriteArray( raster )

outband = outraster.GetRasterBand(1)

#Write to file     
outband.WriteArray(nsidc)
outband.FlushCache()

#Clear arrays and close files
outband = None
iceraster = None
outraster = None
outarray = None

This then is the resulting and imported image:

Friday, February 21, 2014

Reading ERS / ASAR Data and Creating Quicklooks

In this previous post it was described how gdal reads Radarsat SAR data. I want to read ERS /ASAR data and create quicklooks for each file.

ERS /ASAR data has different formats, files ending with *.001 , *.E1, *.E2 or *. N1  Gdal can read all these formats, also the *.001 files which not all software reads. For the latter, the "DAT_01.001" is the file to be passed to gdal.

As first step, I create a list containing all files with these endings. I use os.walk to recursively search through all folders:

filelist_CEOS = []
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, 'DAT_01.001'):
      filelist_CEOS.append(os.path.join(root, filename))

filelist_E1E2 = []
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.E1'):
      filelist_E1E2.append(os.path.join(root, filename))      
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.E2'):
      filelist_E1E2.append(os.path.join(root, filename))  
for root, dirnames, filenames in os.walk('Z:\\ERS_Envisat_SAR\\Arctic\\2005'):
  for filename in fnmatch.filter(filenames, '*.N1'):
      filelist_E1E2.append(os.path.join(root, filename)) 

gdalwarp is able to read all of this formats and I can proceed to process all the files from this filelist, calling gdalwarp from within a Python script. This map-projects the raw SAR file and creates a GeoTIFF:


os.system('gdalwarp -tps  -t_srs EPSG:32633 ' + file_from_filelist + ' ' + outputfilename )  


For a small quicklook a convert it to jpeg and make it smaller:
os.system('gdal_translate -of JPEG -ot byte -outsize 20% 20% -scale 0 1000 0 255 ' + outputfilename + ' ' + browseimage )

The whole script reading in all ERS /ASAR images and creating a quicklook can be found here and is hopefully documented well enough.


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: