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: