Wednesday, October 3, 2012

Creating a Subset of Raster Data with Python / gdal

Continuing to play with Python and gdal...  Reading a GeoTIFF, subsetting it and saving as ENVI file.

First importing the libraries and setting working directory:

   >>> import gdal
   >>> from gdalconst import *
   >>> import os
   >>> os.chdir('C:\Users\max\Documents\PythonProjects')


Opening the file:

   >>> filename = 'ERS1PRI_19920430o04133tr481fr1989_AppOrb_Calib_Spk_SarsimTC_LinDB.tif'
   >>> dataset = gdal.Open(filename, GA_ReadOnly)


I get the band from the opened file. The ReadAsArray with the offset in x and y direction and the x dimensions and y dimensions of the subset:

    >>> band = dataset.GetRasterBand(1)
  >>> subset = band.ReadAsArray(3500, 3500, 50, 50)


The subset can be inspected:

   >>> subset
  array([[-11.49956703, -10.51208401,  -9.58573532, ..., -11.87286568,
        -12.77604961, -15.28282642],
       [-13.6348238 , -14.13398838, -13.61361694, ..., -13.09877014,
        -13.3733263 , -15.26138687],
       [-14.01767826, -12.75454712, -12.39109707, ..., -13.30418301,
        -14.24573517, -15.0796051 ],
       ...,
       [ -9.57748413, -11.63577938, -14.19717503, ...,   3.46974254,
          0.29983696,   0.51451778],
       [-11.11393261, -12.06205177, -13.13371372, ...,   2.4157486 ,
          3.08015943,   0.57344818],
       [-11.68845177, -11.69376755, -11.36115456, ...,   1.95159626,
          1.82419062,   3.39742208]], dtype=float32)
 



Now I register the ENVI driver and Create the ENVI file using the dimensions of the subset:

    >>> driver = gdal.GetDriverByName('ENVI')
  >>> driver.Register()
  >>> outDataset = driver.Create('ERS1PRI_19920430_ENVI_subset', 50, 50, 1, gdal.GDT_Float32)


Finally I get the band of the created ENVI file and write the subset array into the file. After that I can close the files.
   >>> outBand = outDataset.GetRasterBand(1)
   >>> outBand.WriteArray(subset, 0, 0)
   >>> band = None

   >>> dataset = None
   >>> outDataset = None
   >>> subset = None


Closing files is not a "close" command but simply setting the variables to "None"

[I seem to have problems with WriteArray not always writing the file correctly -- unsure yet when and why, the one above works]

No comments:

Post a Comment