Unable to write array to file

development python
1

Hi,
I am trying to use numpy to stack bands and filter the result then write it to file. However I get the error “no product reader” or a “null pointer exception”. I hase seen these documented but no real solution offered. The source file is a stack built with SNAP. I have tried to create a file just with the band, add a band to the existing stack, following the different tutorials but nothing seems to work.
Here is the code.


import os, subprocess
sys.path.append(os.environ['HOME']+os.sep+'.snap'+os.sep+'snap-python')
import esa_snappy as es
import numpy as np

def process_array(p,band, X, Y):
    '''
    Subroutine to extract numpy array from band'''
    print(f"converting band {band} to numpy array")
    source=p.getBand(band)
    width = p.getSceneRasterWidth()
    height = p.getSceneRasterHeight()
    print(f"Array dimensions {height} x {width}")
    arr=np.zeros((height,width))
    source.readPixels(0, 0, width, height, arr)
    arr[arr!=0] -=arr[Y,X]
    print("Array type: ",arr.dtype)
    print("sum of data: ",arr.sum())
    print("Array shape: ", arr.shape)
    return arr

X,Y=2258,1214 #pixel coordinates for referencing
stack=es.ProductIO.readProduct('Path_to_stack.dim')
bands=list(stack.getBandNames())
for i in range(len(bands)):
    arr = process_array(stack, bands[i], X, Y)
    if i ==0:        
        host=np.zeros((len(bands),arr.shape[0],arr.shape[1]))
    host[i] = arr
host_stack=np.sum(host,axis=0)
stack_m=es.Product('Stack_B6_m','Type', stack.getSceneRasterWidth(), stack.getSceneRasterHeight())
es.ProductUtils.copyMetadata(stack, stack_m)
stack_m.setDescription('Calibrated Stack')
stack_m.setProductWriter(es.ProductIO.getProductWriter("BEAM-DIMAP"))
stack_m.writeHeader(es.String('B6_Stack_m.dim'))
es.ProductUtils.copyGeoCoding(stack,stack_m)
targetband=stack_m.addBand('stacked_d', es.ProductData.TYPE_FLOAT32)
targetband.setPixels(0,0,stack.getSceneRasterWidth(),stack.getSceneRasterHeight(),host_stack)
es.ProductIO.writeProduct(stack_m, os.path.join(work_dir,'unwrapped','B6_Stack_m.dim')
                                                    ,"BEAM-DIMAP")

I also tried to
targetband.writePixels(0,0,stack.getSceneRasterWidth(),stack.getSceneRasterHeight(),host_stack)

But it is not happy either. In the examples there are some confusions between row and column oriented arrays. So it is unclear if the array needs to be flat and in which order but I don’t really think that’s the issue, is it?

host_stack is correctly coded and contain data however the layer is only filled with 0s.

Thanks in advance for your answer

2

Have you tried the examples provided by snappy?
Do they work?
Maybe you can use for example the ndvi script as template.
esa-snappy/src/main/resources/esa_snappy/examples/snappy_ndvi.py at master · senbox-org/esa-snappy (github.com)

Yes, the data needs to be flat. Starting from top-left to bottom-right. But as you said this should be the issue in your case.
By looking at your code, I think you can omit the lines
stack_m.setProductWriter(es.ProductIO.getProductWriter("BEAM-DIMAP")) stack_m.writeHeader(es.String('B6_Stack_m.dim'))

But I can’t see obvious errors in your code.

3

Hi,
I have tried all the examples and various methods I saw. A flat array is easy with ravel() at least that’s one thing I know.
What is extremely strange is that it writes the file row by row??? that I have not tried but I cannot see what difference it could make in vectorial computations.
Thanks for the help

4

Also in no example did I see a numpy array object being written directly. This is always with reference with a reader of a band.
Even an array with 0 and 1.

5

In fact if I writePixel before setting the header, I get a pair img/hdr with the array values. I then write the header and it is readable. however when you copy the geocoding, it doesn’t copy the incident angle and slant. So I took it manually from the source file header and copy pasted the relevant information in the header…
This is relatively nuts just to stay polite.

6

Yes, the data of the BEAM-DIMAP format is written line by line. Keep in mind, this format is ~25 years old and hasn’t changed since then. You can switch to NetCDF4-CF or to the latest SNAP format ZNAP. Which is Zarr based.

Incident angle and slant doesn’t really belong to the geocoding.
But if you want to copy most of the data from the source you can use for example ProductUtils.copyProductNodes(), see ProductUtils (SNAP Engine API)

Or one of the specific copy methods:

  • ProductUtils.copyMetadata()
  • ProductUtils.copyTiePointGrids()
  • ProductUtils.copyFlagCodings()
  • ProductUtils.copyFlagBands()
  • ProductUtils.copyGeoCoding()
  • ProductUtils.copyMasks()
  • ProductUtils.copyVectorData()
  • ProductUtils.copyIndexCodings()
  • ProductUtils.copyQuicklookBandName()

For incident and slant range you would use ProductUtils.copyTiePointGrids().

1 Like
7

Hi,
I was referring to the example where the user seems to not know python and write the file row by row. Not sure it is the best to put forward in the python app examples. But since nobody updated the main branch of this repo in 2 years I can imagine that the project is dead?
I managed to write a pair of .img/.hdr files they are just never properly formatted into a BEAM-DIMAP because the header is always poorly formatted by the lib. I edited it by hand and copied the image into the folder to make it work. Not really great…

8

If the writer writes line wise, then it is also good to provide the data line wise. Otherwise, you will experience memory/performance drawbacks.

Probably not dead, but priorities are differently set by ESA at the moment. Might get more attention in the future again, I hope.

The .img/.hdr files serve mainly the BEAM-DIMAP format. They were not intended to be used by users directly in the first place.
But you are right they could be improved. Especially the supported projections. This is also an issue when you use the ENVI writer directly.
What’s your issue with the format specifically?

9

NetCDF4-CF is now well established and supported by many 3rd party tools. Before I retired in 2018 I was using NetCDF4-CF to save products generated outside SNAP and could import them reliably.