DisparityMapToElevationMap failing

OTB Usage
1

I have completed most of the steps to create DEM from a stereo pair of images. The images are 188x1800 pixels at 5 meter resolution (UTM projection). But the final step fails with an error about malloc. And indeed the predicted Elevation map size is HUGE: 40 million x 200 million pixels. I can’t understand where that output size comes from. All preparatory steps completed quickly and produced “normal” sized rasters.

Any tips would be appreciated.

Here’s the command and output:

micha@tp480:GIS$ otbcli_DisparityMapToElevationMap \
>         -io.in ${GIS_dir}BlockMatching_disparity.tif \
>         -io.left ${GIS_dir}evrona_band_5.tif -io.right ${GIS_dir}evrona_band_6.tif \
>         -io.lgrid ${GIS_dir}GridBasedImageResampling_left.tif -io.rgrid ${GIS_dir}GridBasedImageResampling_right.tif \
>         -io.out ${GIS_dir}dem.tif float -elev.dem $SRTM_dir -ram 4196
2019-11-28 00:24:24 (INFO): Default RAM limit for OTB is 128 MB
2019-11-28 00:24:24 (INFO): GDAL maximum cache size is 790 MB
2019-11-28 00:24:24 (INFO): OTB will use at most 8 threads
2019-11-28 00:24:24 (INFO): No kwl metadata found in file ../GIS/BlockMatching_disparity.tif
2019-11-28 00:24:24 (INFO): No kwl metadata found in file ../GIS/evrona_band_5.tif
2019-11-28 00:24:24 (INFO): No kwl metadata found in file ../GIS/evrona_band_6.tif
2019-11-28 00:24:24 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_left.tif
2019-11-28 00:24:24 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_right.tif
2019-11-28 00:24:24 (INFO): Elevation management: setting default height above ellipsoid to 0 meters
2019-11-28 00:24:24 (INFO): Elevation management: using DEM directory (../SRTM/)
2019-11-28 00:24:24 (INFO): Elevation map origin : [691983,3.28782e+06]

2019-11-28 00:24:24 (INFO): Elevation map size : [39733808,200375083]

2019-11-28 00:24:24 (INFO): Estimated memory for full processing: 1.10248e+13MB (avail.: 4196 MB), optimal image partitioning: 2627454039 blocks
2019-11-28 00:24:29 (FATAL): Caught std::exception during application execution: std::bad_alloc

Thanks, Micha

2

I’ve been trying to learn all the required steps for this process myself:

I had similar errors to yours, but only when attempting to use otbcli_StereoFramework. I switched to doing all the individual steps manually as I could not get that working.

I notice you don’t have kwl metadata found in your files.

Are your source files from one of these supported sensors:

  • Pleiades
  • SPOT5
  • QuickBird
  • Ikonos
  • WorldView-1
  • WorldView-2
  • Formosat

If they are you should have that data in a linked file that came with the source data.

I found that things did not go smoothly, unless I maintained the kwl metadata against my rasters. If I somehow lost that metadata through the various processes, I went back and re-worked my steps.

I certainly have not nailed the process down yet as you can see from my results, it would be great to see some screen shots of your output?

3

Thanks, Rossoe, for your reply

I am using bands from the Venus satellite.
https://venus.bgu.ac.il/venus/
I guess there is no keyword list file with this data. Any idea where I can find the structure of the kwl metadata file?

4

Here’s one of the RPB’s from one of my (NEXTVIEW) stereo pair rasters, to give you an idea of the construct:
16JUN18143537-P1BS-057461431010_01_P001.RPB (1.8 KB)

I’m not sure how easy it is to generate these yourself.

I think it’s either a RPB or RPC file that you would be looking for.

This article describes the variations in type of RPC files you get with different with different sensor models - https://www.harrisgeospatial.com/docs/RPCOrthorectification.html

5

Thanks for that help.
I found that the satellite data I’m using has *.HDR files. So when I cropped out my ROI (using ExtractROI) I was careful to give the cropped raster the same name as the original, and I copied the HDR to the crop subdirectory. Then in all subsequent processes, the “No kwl file” warnings disappeared.

However, when I get to the final step, DisparityMapToElevationMap, I still see that extremely large predicted DSM size. And the process fails since it cannot allocate such large memory.
I have no idea how to proceed from here :anguished:

6

Could you tell us how you are cropping your ROI?

I found that using the ‘standard’ or ‘extent’ options just didn’t work for me. So I switched to using fit (vector) and fed in a shapefile of my desired ROI.

-mode.fit.vect vectorfile

So my whole ExtractROI command was:

otbcli_ExtractROI -in F:/GIS/02AUG12WV020800012AUG02114451-P1BS-052766399070_01_P001.ntf -out F:/GIS/14_11_19/stereo_cliff_2_small_Ls.tif -mode fit -mode.fit.vect F:/GIS/square_final.shp -ram 128

Also when you finish the otbcli_StereoRectificationGridGenerator step what dimensions do you get reported at the end?

Should look something like:

Output parameters value:
epi.rectsizex: 5908
epi.rectsizey: 5907
epi.baseline: 0.3019381762

Are you using command line, or the GUI interface to complete these steps?

7

Here’s my extract command:

otbcli_ExtractROI -in ${VENUS_file} -out ${evrona_file} -mode standard -startx 5030 -starty 1220 -sizex 800 -sizey 1800 2019-11-29 15:20:33 (INFO): No kwl metadata found in file ../VE_VM01_VSC_L1VALD_ISRAES10_20190306.DBL.DIR/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306.DBL.TIF
2019-11-29 15:20:33 (INFO): Default RAM limit for OTB is 128 MB
2019-11-29 15:20:33 (INFO): GDAL maximum cache size is 790 MB
2019-11-29 15:20:33 (INFO): OTB will use at most 8 threads
2019-11-29 15:20:33 (INFO): Estimated memory for full processing: 329.018MB (avail.: 128 MB), optimal image partitioning: 3 blocks
2019-11-29 15:20:33 (INFO): File ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306.DBL.TIF will be written in 4 blocks of 800x451 pixels
Writing ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306.DBL.TIF...: 100% [**************************************************] (0 seconds)

I used -mode standard by setting the row-column pixel values (not real projected coordinates)

And here’s the output of StereoRectification:

micha@tp480:code$ otbcli_StereoRectificationGridGenerator -io.inleft ${band_5} -io.inright ${band_6} -io.outleft ${GIS_dir}StereoRectificationGridGenerator_left.tif float -io.outright ${GIS_dir}StereoRectificationGridGenerator_right float -epi.elevation.dem ${SRTM_dir} -epi.elevation.default 10
2019-11-29 15:32:48 (INFO): Default RAM limit for OTB is 128 MB
2019-11-29 15:32:48 (INFO): GDAL maximum cache size is 790 MB
2019-11-29 15:32:48 (INFO): OTB will use at most 8 threads
2019-11-29 15:32:48 (INFO): No kwl metadata found in file ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306_5.DBL.TIF
2019-11-29 15:32:48 (INFO): No kwl metadata found in file ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306_6.DBL.TIF
2019-11-29 15:32:48 (INFO): Elevation management: setting default height above ellipsoid to 10 meters
2019-11-29 15:32:48 (INFO): Elevation management: using DEM directory (../SRTM/)
Computing epipolar grids ...: 100% [**************************************************] (0 seconds)
2019-11-29 15:32:48 (INFO): Estimated memory for full processing: 1626.49MB (avail.: 128 MB), optimal image partitioning: 13 blocks
2019-11-29 15:32:48 (INFO): File ../GIS/StereoRectificationGridGenerator_left.tif will be written in 18 blocks of 336x336 pixels
Writing ../GIS/StereoRectificationGridGenerator_left.tif...: 4% [**                                                ]2019-11-29 15:32:48 (WARNING): Skipping GCPs saving to prevent GDAL from assigning a WGS84 projref to file (../GIS/StereoRectificationGridGenerator_left.tif)
Writing ../GIS/StereoRectificationGridGenerator_left.tif...: 100% [**************************************************] (0 seconds)
2019-11-29 15:32:48 (WARNING): Check filename: no extension detected, using TIF as default.
2019-11-29 15:32:48 (INFO): Estimated memory for full processing: 50.4891MB (avail.: 128 MB), optimal image partitioning: 1 blocks
2019-11-29 15:32:48 (INFO): File ../GIS/StereoRectificationGridGenerator_right.tif will be written in 1 blocks of 1802x802 pixels
Writing ../GIS/StereoRectificationGridGenerator_right.tif...: 100% [**************************************************]2019-11-29 15:32:48 (WARNING): Skipping GCPs saving to prevent GDAL from assigning a WGS84 projref to file (../GIS/StereoRectificationGridGenerator_right.tif)
 (0 seconds)
Output parameters value:
epi.rectsizex: 1800
epi.rectsizey: 800
epi.baseline: 0

And I see that the “no kwl metadata” warnings are still persisting…

8

I was confused for a minute, as you mentioned dimensions of 188x1800 in first post, but in ExtractROI you are cropping to 1800 x 800 dimensions.

What does the output raster look like at the otbcli_DisparityMapToElevationMap step?

9

Yes, I had a typo error. My ROI is 800x1800.
No output raster from DisparityMapToElevation:

micha@tp480:code$ otbcli_DisparityMapToElevationMap \
>         -io.in ${GIS_dir}BlockMatching_disparity.tif \
>         -io.left ${band_5} -io.right ${band_6} \
>         -io.lgrid ${GIS_dir}GridBasedImageResampling_left.tif -io.rgrid ${GIS_dir}GridBasedImageResampling_right.tif \
>         -io.out ${Out_dir}evrona_dsm.tif float -elev.dem $SRTM_dir -ram 4196
2019-12-01 08:32:20 (INFO): Default RAM limit for OTB is 128 MB
2019-12-01 08:32:20 (INFO): GDAL maximum cache size is 790 MB
2019-12-01 08:32:20 (INFO): OTB will use at most 8 threads
2019-12-01 08:32:20 (INFO): No kwl metadata found in file ../GIS/BlockMatching_disparity.tif
2019-12-01 08:32:20 (INFO): No kwl metadata found in file ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306_5.DBL.TIF
2019-12-01 08:32:20 (INFO): No kwl metadata found in file ../GIS/VE_VM01_VSC_PDTIMG_L1VALD_ISRAES10_20190306_6.DBL.TIF
2019-12-01 08:32:20 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_left.tif
2019-12-01 08:32:20 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_right.tif
2019-12-01 08:32:20 (INFO): Elevation management: setting default height above ellipsoid to 0 meters
2019-12-01 08:32:20 (INFO): Elevation management: using DEM directory (../SRTM/)
2019-12-01 08:32:20 (INFO): Elevation map origin : [691983,3.28782e+06]

2019-12-01 08:32:20 (INFO): Elevation map size : [39733808,200375083]

2019-12-01 08:32:20 (INFO): Estimated memory for full processing: 1.10248e+13MB (avail.: 4196 MB), optimal image partitioning: 2627454039 blocks
2019-12-01 08:32:30 (FATAL): Caught std::exception during application execution: std::bad_alloc

Note the line: Elevation map size: I can’t figure out where those huge sizes are coming from.

Here’s another weird result. After running the otbcli_StereoRectificationGridGenerator, I get the two rectification grids with this info:

micha@tp480:code$ otbcli_ReadImageInfo -in ../GIS/StereoRectificationGridGenerator_left.tif 
2019-12-01 08:40:18 (INFO): Default RAM limit for OTB is 128 MB
2019-12-01 08:40:18 (INFO): GDAL maximum cache size is 790 MB
2019-12-01 08:40:18 (INFO): OTB will use at most 8 threads
2019-12-01 08:40:18 (INFO): No kwl metadata found in file ../GIS/StereoRectificationGridGenerator_left.tif
2019-12-01 08:40:18 (INFO): 
Image general information:
        Number of bands : 2
        No data flags : Not found
        Start index :  [0,0]
        Size :  [1802,802]
        Origin :  [0,0]
        Spacing :  [5,5]
        Estimated ground spacing (in meters): [7147.66,28881.8]

Image acquisition information:
        Sensor : 
        Image identification number: 

Image default RGB composition:
        [R, G, B] = [0,1,2]

Ground control points information:
        Number of GCPs = 0
        GCPs projection = 

Output parameters value:
indexx: 0
indexy: 0
sizex: 1802
sizey: 802
spacingx: 5
spacingy: 5
originx: 0
originy: 0
estimatedgroundspacingx: 7147.65625
estimatedgroundspacingy: 28881.78906
numberbands: 2
sensor: 
id: 
time: 
town: 
country: 
rgb.r: 0
rgb.g: 1
rgb.b: 2
projectionref: 
keyword: 
gcp.count: 0
gcp.proj: 
gcp.ids: 
gcp.info: 
gcp.imcoord: 
gcp.geocoord:

The size of the rectification grids looks OK [1802, 802], but the ground spacing values are totally off.
Shouldn’t the spacing be about 4000x9000 (five times the grid size)???

Regards, Micha

10

Sorry my typo, I meant to say what dimensions are you getting at the otbcli_BlockMatching step?

Interested to know if dimensions are correct right up to BlockMatching step, but then only go rogue on the final otbcli_DisparityMapToElevationMap step?

I’m not sure about the expectations re ground spacing, that’s beyond my knowledge unfortunately - hopefully someone else on here might know.

For reference here is the results for me using otbcli_ReadImageInfo on one of my rectification grids:

F:\OTB-6.6.1-Win64\bin>otbcli_ReadImageInfo -in F:/GIS/14_11_19/31_grid_stereo_2_small_L2.tif
2019-12-02 02:11:35 (INFO): Default RAM limit for OTB is 128 MB
2019-12-02 02:11:35 (INFO): GDAL maximum cache size is 812 MB
2019-12-02 02:11:35 (INFO): OTB will use at most 8 threads
2019-12-02 02:11:35 (INFO): No kwl metadata found in file F:/GIS/14_11_19/31_grid_stereo_2_small_L2.tif
2019-12-02 02:11:35 (INFO):
Image general information:
        Number of bands : 2
        No data flags : Not found
        Start index :  [0,0]
        Size :  [5760,5756]
        Origin :  [0,0]
        Spacing :  [1,1]
        Estimated ground spacing (in meters): [2984.31,2417.28]

Image acquisition information:
        Sensor :
        Image identification number:

Image default RGB composition:
        [R, G, B] = [0,1,2]

Ground control points information:
        Number of GCPs = 0
        GCPs projection =

Output parameters value:
indexx: 0
indexy: 0
sizex: 5760
sizey: 5756
spacingx: 1
spacingy: 1
originx: 0
originy: 0
estimatedgroundspacingx: 2984.313477
estimatedgroundspacingy: 2417.278809
numberbands: 2
sensor:
id:
time:
town:
country:
rgb.r: 0
rgb.g: 1
rgb.b: 2
projectionref:
keyword:
gcp.count: 0
gcp.proj:
gcp.ids:
gcp.info:
gcp.imcoord:

gcp.geocoord:
11

Hello Ross:

Here are the command and image info of my BlockMatching step. Does anything “jump out” at your??
(Note that the BlockMatching output is set to 300X600 pixels. ??
and the origin is (0,0) even tho’ the CRS is still UTM 36N ??)

micha@tp480:code$ otbcli_BlockMatching \
>         -io.inleft ${GIS_dir}GridBasedImageResampling_left.tif \
>         -io.inright ${GIS_dir}GridBasedImageResampling_right.tif \
>         -io.out ${GIS_dir}BlockMatching_disparity.tif float \
>         -bm.minhd 45 -bm.maxhd 55 -bm.minvd 1 -bm.maxvd -1 io.outmetric 1 \
>         -bm.subpixel triangular -ram 4196
2019-12-08 20:39:43 (INFO): Default RAM limit for OTB is 128 MB
2019-12-08 20:39:43 (INFO): GDAL maximum cache size is 790 MB
2019-12-08 20:39:43 (INFO): OTB will use at most 8 threads
2019-12-08 20:39:43 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_left.tif
2019-12-08 20:39:43 (INFO): No kwl metadata found in file ../GIS/GridBasedImageResampling_right.tif
2019-12-08 20:39:43 (INFO): Estimated memory for full processing: 9.76334MB (avail.: 4196 MB), optimal image partitioning: 1 blocks
2019-12-08 20:39:43 (INFO): File ../GIS/BlockMatching_disparity.tif will be written in 1 blocks of 300x600 pixels
SSD block matching: 100% [**************************************************] (0 seconds)            ]
Sub-pixel refinement: 100% [**************************************************] (0 seconds)
Writing ../GIS/BlockMatching_disparity.tif...: 100% [**************************************************] (0 seconds)


micha@tp480:code$ otbcli_ReadImageInfo -in ../GIS/BlockMatching_disparity.tif
2019-12-08 20:41:39 (INFO): Default RAM limit for OTB is 128 MB
2019-12-08 20:41:39 (INFO): GDAL maximum cache size is 790 MB
2019-12-08 20:41:39 (INFO): OTB will use at most 8 threads
2019-12-08 20:41:39 (INFO): No kwl metadata found in file ../GIS/BlockMatching_disparity.tif
2019-12-08 20:41:39 (INFO): 
Image general information:
        Number of bands : 2
        No data flags :
                Band 1: 0
                Band 2: 0
        Start index :  [0,0]
        Size :  [300,600]
        Origin :  [0,0]
        Spacing :  [5,5]
        Estimated ground spacing (in meters): [4.9804,5.01458]

Image acquisition information:
        Sensor : 
        Image identification number: 
        Image projection : PROJCS["WGS 84 / UTM zone 36N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",33],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32636"]]

Image default RGB composition:
        [R, G, B] = [0,1,2]

Ground control points information:
        Number of GCPs = 0
        GCPs projection = 

Output parameters value:
indexx: 0
indexy: 0
sizex: 300
sizey: 600
spacingx: 5
spacingy: 5
originx: 0
originy: 0
estimatedgroundspacingx: 4.980398655
estimatedgroundspacingy: 5.014582157
numberbands: 2
sensor: 
id: 
time: 
town: 
country: 
rgb.r: 0
rgb.g: 1
rgb.b: 2
projectionref: PROJCS["WGS 84 / UTM zone 36N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",33],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32636"]]
keyword: 
gcp.count: 0
gcp.proj: 
gcp.ids: 
gcp.info: 
gcp.imcoord: 
gcp.geocoord:

Thanks for following up on this.

12

Just had a thought on this one, you mentioned :

when I get to the final step, DisparityMapToElevationMap, I still see that extremely large predicted DSM size. And the process fails since it cannot allocate such large memory.
I have no idea how to proceed from here

Have you tried entering a higher -step value into the otbcli_DisparityMapToElevationMap command?

If you try -step 2 or -step 3 you might get around the memory issue, as it reduces the resolution of the output DEM.