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:

No comments:

Post a Comment