[gdal-dev] JP2k, ReadError, Python

[Even Rouault]

Le samedi 30 juillet 2011 18:46:42, Jay L. a écrit :
> 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["fal
> se_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.

Ah ok, so it is the Jasper JPEG2000 driver that was used. Well it is known not 
to be very robust, so I'm not surprised that it fails decoding your image.

> 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.

Hum, well then this image is quite particular then if it also fails with the 
JP2ECW driver...

> 
> 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.

Yes, upload it somewhere and provide the link.

> 
> 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?

ok, I understand better. GetMinimum() and GetMaximum() do zero computation. 
They only return a value if it is already known. If it is not known, then 
ComputeRasterMinMax() is definitely the way to go.

> 
> 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