I am having an issue reading an ENVI binary file (with .hdr file) from a zip file using GDAL in Python and the /vsizip driver.
I know the ENVI file (binary and header) itself is fine, since when I unzip the binary file and the .hdr file from the archive and use gdal_translate directly it works fine.
gdal_translate phsig.cor.geo test.tif
Input file size is 5139, 6060
0...10...20...30...40...50...60...70...80...90...100 - done.
If I use gdal_translate to attempt the same operation with the zip file I get the following result.
gdal_translate /vsizip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf.zip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf/merged/phsig.cor.geo test2.tif
Input file size is 3681, 4083
0...10...20...30...40...50...60...70...80...90...100 - done.
The dimensions reported above (3681,4083) are wrong. The file is actually (5139, 6060). I have not been able to find a good discussion of the error shown above.
I also have implemented effectively this same function in python. When I open and read the file from my Python script I am able to get some information correctly, such as number of bands. However, the image dimensions are not read correctly and so the data are broadcast into an array that is the wrong dimension.
My function for reading is:
def __init__(self,imName):
handle = gdal.Open(imName)
self.numBands = handle.RasterCount
self.geoTransform = handle.GetGeoTransform()
self.projection = handle.GetProjection()
self.x = handle.RasterXSize
self.y = handle.RasterYSize
self.data = np.zeros((self.y,self.x,self.numBands))
self.min = np.zeros(self.numBands)
self.max = np.zeros(self.numBands)
self.bandType = []
for i in range(0,self.numBands):
banddata = handle.GetRasterBand(i+1)
self.data[:,:,i] = banddata.ReadAsArray()
self.bandType.append(gdal.GetDataTypeName(banddata.DataType).lower())
self.min[i] = np.min(self.data[:,:,i])
self.max[i] = np.max(self.data[:,:,i])
if self.min[i] is None or self.max[i] is None:
(self.min[i],self.max[i]) = banddata.ComputeRasterMinMax(1)
self.metadata = handle.GetMetadata()
self.gcp = handle.GetGCPs()
self.gcpProj = handle.GetGCPProjection()
It appears that gdal is not reading the .hdr file properly from the zip archive. However, it seems this functionality should be available from looking at the gdal documentation and the specification for the vsizip driver.
Following up from Luke's comments, my GDAL version is 2.1.2, verified with gdalinfo --version. The output from gdalinfo for the file in unzipped form and within the zip archive are below.
/opt/local/bin/gdalinfo -stats phsig.cor.geo
Driver: ENVI/ENVI .hdr Labelled
Files: phsig.cor.geo
phsig.cor.geo.hdr
Size is 5139, 6060
Coordinate System is:
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-112.731111111111105,33.272777777777776)
Pixel Size = (0.000277777777778,-0.000277777777778)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-112.7311111, 33.2727778) (112d43'52.00"W, 33d16'22.00"N)
Lower Left (-112.7311111, 31.5894444) (112d43'52.00"W, 31d35'22.00"N)
Upper Right (-111.3036111, 33.2727778) (111d18'13.00"W, 33d16'22.00"N)
Lower Right (-111.3036111, 31.5894444) (111d18'13.00"W, 31d35'22.00"N)
Center (-112.0173611, 32.4311111) (112d 1' 2.50"W, 32d25'52.00"N)
Band 1 Block=5139x1 Type=Float32, ColorInterp=Undefined
Minimum=0.000, Maximum=1.000, Mean=0.591, StdDev=0.476
Metadata:
STATISTICS_MAXIMUM=1
STATISTICS_MEAN=0.590579634219
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=0.47649836193737
/opt/local/bin/gdalinfo -stats /vsizip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf.zip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf/merged/phsig.cor.geo
Driver: ENVI/ENVI .hdr Labelled
Files: /vsizip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-
76d1-v1.1.3-asf.zip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf/merged/phsig.cor.geo
/vsizip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf.zip/S1-IFG_RM_M1S2_TN027_20170301T132816-20170217T132724_s2-poeorb-76d1-v1.1.3-asf/merged/phsig.cor.hdr
Size is 3681, 4083
Coordinate System is:
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0000000, 0.0000000) ( 0d 0' 0.01"E, 0d 0' 0.01"N)
Lower Left ( 0.000, 4083.000) ( 0d 0' 0.01"E,Invalid angle)
Upper Right ( 3681.000, 0.000) (Invalid angle, 0d 0' 0.01"N)
Lower Right ( 3681.000, 4083.000) (Invalid angle,Invalid angle)
Center ( 1840.500, 2041.500) (Invalid angle,Invalid angle)
Band 1 Block=3681x1 Type=Float32, ColorInterp=Undefined
Minimum=0.000, Maximum=1.000, Mean=0.575, StdDev=0.471
Metadata:
STATISTICS_MAXIMUM=1
STATISTICS_MEAN=0.575013162695
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=0.4709390447683
From this it appears that while reading the file using /vsizip it is not reading the geotranform block either. The projection information seems correct, but the image coordinates are not converted to geographic coordinates properly. This is consistent with the behavior from Python as well.
