[gdal-dev] JP2k, ReadError, Python

[Jay L.]

Evan,

Thanks for the info.  Setting debug to on has allowed me to trace the
problem to this line:

scanBlock = band.ReadAsArray(h,i,numberColumns,
numberRows).astype(numpy.float)

showing this error / information:

C:\Py_Scripts\RunningSD>python ISCopy.py P01_001468_1534_XI_26S034W.jp2 -vus
5
JPEG2000: IHDR box found. Dump: width=9973, height=34577, numcmpts=1, bpp=8
JPEG2000: Component 0: bpp=8, signedness=0
GDALJP2Metadata: Got projection from GeoJP2 (geotiff) box:
PROJCS["Mars_Equidistant_Cylindrical_Sphere",GEOGCS["GCS_Mars_Sphere",DATUM["Mars_Sphere",SPHEROID["Mars_Spher
METER["standard_parallel_1",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
GDAL: GDALOpen(P01_001468_1534_XI_26S034W.jp2, this=034E8230) succeeds as
JPEG2000.
GDAL: QuietDelete(output.tif) invoking Delete()
GDAL: GDALOpen(output.tif, this=01D5F260) succeeds as GTiff.
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALClose(output.tif, this=01D5F260)
GDAL: GDALDriver::Create(GTiff,output.tif,9973,34577,1,Byte,00000000)
Processing band: 1
error: cannot decode code stream
JPEG2000: Unable to decode image. Format: jp2, 4
Traceback (most recent call last):
  File "ISCopy.py", line 357, in <module>
    main()
  File "ISCopy.py", line 343, in main
    perform_stretch(dataset,outdataset,**options.__dict__)
  File "ISCopy.py", line 235, in perform_stretch
    input_array = normalize(input_array, bandMin, bandMax)
  File "ISCopy.py", line 52, in normalize
    inputBand -= bandMin
  File "C:\Python26\ArcGIS10.0\lib\site-packages\numpy\ma\core.py", line
3698, in __isub__
    ndarray.__isub__(self._data, np.where(self._mask, 0, getdata(other)))
TypeError: unsupported operand type(s) for -: 'NoneType' and 'float'
GDAL: GDALClose(P01_001468_1534_XI_26S034W.jp2, this=034E8230)
GDAL: GDALClose(output.tif, this=01D5F260)
GDAL: GDALDeregister_GTiff() called.

Removing "astype" solves the attribute error, but the line is still
returning "None" when it should be returning an array.  Entering
gdal_translate in the command line without arguments shows that:

The following format drivers are configured and support output:
  ECW: ERDAS Compressed Wavelets (SDK 3.x)
  JP2ECW: ERDAS JPEG2000 (SDK 3.x)

Using "gdal.GetDriverByName("JPEG2000").Deregister()" to unload the JPEG2000
driver results in this error:

RuntimeError: `P01_001468_1534_XI_26S034W.jp2' not recognised as a supported
file format.

The ECW driver was downloaded and installed from gisinternals last evening.

I have a sample JP2 that is 28MB zipped.  What is the best way to get that
to you?  I can upload it to a webspace and get you the link if that is
easiest.

Not sure to follow you. What is the  "fast" method for getting the min/max ?
> band.ComputeRasterMinMax(False) ?


I had found that band.GetMinimum() and band.GetMaximum(), which processed
quite quickly, often returned "None" is statistics had not already been
calculated for the raster.  Using band.ComputeRasterMinMax(), while slower,
always allowed me to get a value.  Is there a better way to get the band
Min/Max?

Many thanks,
Jay


On Sat, Jul 30, 2011 at 9:12 AM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> Le samedi 30 juillet 2011 17:37:48, Jay L. a écrit :
> > Morning all,
> >
> > I have a rather large JP2 file (created using the Kakadu driver I
> believe)
> > that I am trying to process with a python script.  I am able to open the
> > dataset and grab the first band.
> >
> > To calculate band statistics I use the following code which returns 1000
> > lines of error before exiting
> >
> > def get_band_stats(band):
> >     NDV = band.GetNoDataValue()
> >     bandMin = band.GetMinimum()
> >     bandMax = band.GetMaximum()
> >     if bandMin is None or bandMax is None:
> >         #Approx has to be set to false to ensure that an accurate Min/Max
> > are calculated...
> >         (bandMin, bandMax) = band.ComputeRasterMinMax(False)
> >     return NDV, bandMin, bandMax
> >
> > ....
> > ERROR 1: IReadBlock failed at X offset 20, Y offset 25
> > ERROR 1: IReadBlock failed at X offset 21, Y offset 25
> > ....
> >
> > This cuts off after 1000 lines.
> >
> > Commenting out the if statement, thereby using the "fast" method for
> > getting the min/max returns the following error:
> >
> > error: cannot decode code stream
> > ERROR 1: IReadBlock failed at X offset 0, Y offset 0
> > ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0
>
> Not sure to follow you. What is the  "fast" method for getting the min/max
> ?
> band.ComputeRasterMinMax(False) ?
>
> Did you try running with CPL_DEBUG=ON defined as an environment variable to
> check which driver is used and if there are other interesting warnings ?
>
> >
> > Finally, using gdalinfo I see that the image is being opened with the
> > ECWJP2000 driver.  I am also able to calculate statistics without a
> problem
> > and the returned values are within the expected range.
>
> Hum, strange, ComputeRasterMinMax() and ComputeRasterStatistics() are very
> silmilar : they use the same way of fetching the data, so I'd assume that
> if
> one works, the other one would too.
>
> > I am also able to
> > use gdal_translate to convert these to .tif prior to processing, but I
> > would rather avoid this method.  Do I need to alter my code to work with
> > JP2 as they are compressed?
>
> Theroretically no. If it doesn't work, it's a defect of the driver. Could
> you
> share your JP2 file ?
>
> >
> > Can I implicitly define the driver to be used with JP2 within my python
> > code without performing the unload all / load what I want method?
>
> You can unload a driver this way :
> gdal.GetDriverByName("short_name_of_driver").Deregister()
>
> >
> > Many thanks,
> > Jay
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20110730/a3ac27de/attachment-0001.html