2

I am trying to read an ENVI file as an array using GDAL and Python

Image info are following:

Driver: ENVI/ENVI .hdr Labelled
Files: IMG-VV-ALPSRP248213250-P1.1__D_pwr_geo_sigma
     IMG-VV-ALPSRP248213250-P1.1__D_pwr_geo_sigma.hdr
Size is 1659, 2775
Coordinate System is:
PROJCS["UTM Zone 16, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-87],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (332125.000000000000000,2017650.000000000000000)
Pixel Size = (25.000000000000000,-25.000000000000000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  332125.000, 2017650.000) ( 88d35'16.06"W, 18d14'29.96"N)
Lower Left  (  332125.000, 1948275.000) ( 88d34'55.98"W, 17d36'53.41"N)
Upper Right (  373600.000, 2017650.000) ( 88d11'44.12"W, 18d14'40.22"N)
Lower Right (  373600.000, 1948275.000) ( 88d11'29.00"W, 17d37' 3.30"N)
Center      (  352862.500, 1982962.500) ( 88d23'21.23"W, 17d55'47.10"N)
Band 1 Block=1659x1 Type=Float32, ColorInterp=Undefined

My code is the following:

driver = gdal.GetDriverByName('ENVI')
driver.Register()

#Mind the suffix (It is an ENVI file)    
file = 'C:/img.1__A_pwr_geo_sigma

raster = gdal.Open(file,gdal.GA_ReadOnly)

raster_array = raster.ReadAsArray()

print raster_array

Output:

>>>array([[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

All values within the array are NaN although I know that there are floating point 32 bit values within the image (Checked with ENVI software)

What am I doing wrong here? Or is there a problem with the suffix?

I have also tried translating the ENVI format to Geotiff using gdal_translate but also the GeoTiff produces the same array.

Matthew
  • 9,851
  • 4
  • 46
  • 77
Dimitris
  • 485
  • 1
  • 4
  • 16
  • Could be a limitation with the ENVI driver. What version of GDAL? Is there anywhere the raster can be downloaded? – Mike T Feb 01 '13 at 19:50
  • Thanks for your reply. GDAL version is 1.9.2 Raster is huge, I will clip part of it and upload it to a web location if no other ideas about the problem come up – Dimitris Feb 01 '13 at 21:39
  • Mike can you still help? I have prepared a link in Google drive for you. How can I send you a private message with the link? – Dimitris Feb 02 '13 at 15:17
  • InMail feature of LinkedIn should work; see url on my profile here to get there. – Mike T Feb 02 '13 at 17:57

1 Answers1

2

Good news, not all the values are nan. The GDAL driver (v.1.9.1) works great, and valid data is visible in software that depends on GDAL (e.g., QGIS).

The nan values are used to represent NoData. (Normally, most folks use some finite "dummy" value for NoData, e.g. 1e+20).

import numpy as np
from osgeo import gdal
fname = r'C:\path\to\envi_data\IMG-HH-ALPSRP136260330-H1.1__A_pwr_geo_sigma'
ds = gdal.Open(fname, gdal.GA_ReadOnly)
band = ds.GetRasterBand(1)
print band.GetNoDataValue()
# None ## normally this would have a finite value, e.g. 1e+20
ar = band.ReadAsArray()
print np.isnan(ar).all()
# False
print '%.1f%% masked' % (np.isnan(ar).sum() * 100.0 / ar.size)
# 43.0% masked

The best way to represent arrays with missing values in Python/NumPy is to use a masked array:

mar = np.ma.masked_array(ar, np.isnan(ar))
print mar.min(), np.median(mar), mar.mean(), mar.std(), mar.max()
# 0.000715672 0.148093774915 0.0921848740388 0.0700106260235 5.0
Mike T
  • 41,085
  • 18
  • 152
  • 203