[gdal-dev] gdalwarp with RPCs causing bogus output

Mon Apr 22 07:30:40 PDT 2024

Hi Joe,

Now you've pointed them out, I do see those little bits in the north east.
I hadn't zoomed in enough on the right area. Sadly, I've not much advice to
offer on whether it's expected, and/or how to avoid it.

(Apologies for the slow and unhelpful reply!)

Cheers,
Daniel

On Thu, 18 Apr 2024, 19:02 Joseph McGlinchy, <jmcglinchy at hydrosat.com>
wrote:

> Hi Daniel,
>
> Thank you.  I tried your commands exactly with gdal 3.8.3 and am seeing
> the odd type of output I was describing in the northeastern most reaches of
> the images. Would you be willing to verify if you see that or not?
>
> Here is the result of the debug messaging. I think a screenshot which
> shows what I see results in the message being over the size limit.
>
> GDAL: GDALOpen(../manyPts4000_SRTM_test_NptRPC.tiff, this=0x5584b5c21370)
> succeeds as GTiff.
> GDAL: GDALOpen(output_gdal38.tiff, this=0x5584b5c4c950) succeeds as GTiff.
> GDALWARP: Defining SKIP_NOSOURCE=YES
> Processing ../manyPts4000_SRTM_test_NptRPC.tiff [1/1] : 0WARP: Copying
> metadata from first source to destination dataset
> WARP: Set SOURCE_EXTRA=5 warping options due to RPC warping
> WARP: Set SAMPLE_STEPS=ALL warping options due to RPC DEM warping
> GDAL: GDALOpen(srtm_12_05.tif, this=0x5584b5c522a0) succeeds as GTiff.
> RPC: Short-circuiting coordinate transformation from DEM SRS to WGS 84 due
> to apparent nop
> GDAL: GDAL_CACHEMAX = 1504 MB
> GTiff: ScanDirectories()
> GDAL: GDALDefaultOverviews::OverviewScan()
> Using internal nodata values (e.g. 0) for image
> ../manyPts4000_SRTM_test_NptRPC.tiff.
> WARP: band=0 dstNoData=0.000000
> RPC: Failed Iterations 101: Got: -121.6743548438707,37.81305392714289
>  Offset=0.200509,-7.35862
> WARP: At least one point failed after direct transform
> RPC: Failed Iterations 101: Got: -121.6743548438707,37.81305392714289
>  Offset=0.200509,-7.35862
> RPC: Failed Iterations 101: Got: -121.6756957607436,37.81515846637773
>  Offset=-0.0317686,-12.1573
> RPC: Failed Iterations 101: Got: -121.6759688314497,37.81397923057818
>  Offset=-0.0782949,-21.942
> RPC: Failed Iterations 101: Got: -121.6760567085541,37.81393132082147
>  Offset=-0.0200723,-13.2144
> RPC: Failed Iterations 101: Got: -121.6758750148136,37.8150256294207
>  Offset=0.0121725,-8.26664
> RPC: Failed Iterations 101: Got: -121.6758055299304,37.81513471897441
>  Offset=0.0111142,-9.29081
> RPC: Failed Iterations 101: Got: -121.6757471773221,37.81517816256341
>  Offset=0.00786933,-9.66179
> RPC: Failed Iterations 101: Got: -121.6757610763044,37.81516845364766
>  Offset=0.00852169,-9.58453
> GDAL: GDALSuggestedWarpOutput(): 1 out of 3249 points failed to transform.
> GDAL: GDALWarpKernel()::GWKNearestShort() Src=0,0,2383x2800
> Dst=0,0,2411x3318
> ...10...20...30...40...50GDAL: GDALWarpKernel()::GWKNearestShort()
> Src=1710,0,1460x1401 Dst=2411,0,1205x1659
> ...60GDAL: GDALWarpKernel()::GWKNearestShort() Src=2833,0,1263x2458
> Dst=3616,0,1206x829
> ...GDAL: GDALWarpKernel()::GWKNearestShort() Src=2995,148,1101x1054
> Dst=3616,829,1206x830
> 70..GDAL: GDALWarpKernel()::GWKNearestShort() Src=2036,993,2060x1807
> Dst=2411,1659,2411x1659
> .80...90...100 - done.
> GDAL: GDALClose(srtm_12_05.tif, this=0x5584b5c522a0)
> GDAL: Flushing dirty blocks:
> 0...10...20...30...40...50...60...70...80...90...100 - done.
> GDAL: GDALClose(output_gdal38.tiff, this=0x5584b5c4c950)
> GDAL: GDALClose(../manyPts4000_SRTM_test_NptRPC.tiff, this=0x5584b5c21370)
> ------------------------------
> *From:* Daniel Evans <daniel.fred.evans at gmail.com>
> *Sent:* Thursday, April 18, 2024 9:03 AM
> *To:* Joseph McGlinchy <jmcglinchy at hydrosat.com>
> *Cc:* Even Rouault <even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
> <gdal-dev at lists.osgeo.org>
> *Subject:* Re: [gdal-dev] gdalwarp with RPCs causing bogus output
>
> Hi Joe,
>
> My minimal change is just to add the missing "RPC_DEM=srtm_12_05.tif" and
> remove the "RPC_HEIGHT=350". RPC_HEIGHT appears to override the use of any
> supplied DEM, and the docs at least hint that this is the expected
> behaviour.
>
> ```
> gdalwarp --config CPL_DEBUG ON \
>     -t_srs EPSG:32610 \
>     -rpc \
>     -to "RPC_DEM_SRS=+proj=longlat +datum=WGS84 +no_def" \
>     -to "RPC_DEM_MISSING_VALUE=0.001" \
>     -to "RPC_FOOTPRINT='POLYGON ((-123.12941327726458 38.01872167010621,
> -121.67068342817124 37.82615456858018, -121.87575337280391
> 37.11164114368971, -123.32135086275996 37.30124603928108,
> -123.12941327726458 38.01872167010621))'" \
>     -to "RPC_MAX_ITERATIONS=101" \
>     -to "RPC_DEM=srtm_12_05.tif" \
>     manyPts4000_SRTM_test_NptRPC.tiff \
>     output.tiff
> ```
>
> Changing to "RPC_DEM_SRS=EPSG:4326+5773" as Even suggested gives a
> noticeable change around the DEM locations of dams and lakes with the
> updated CRS, e.g. west and north of Livermore CA. Whether that's an
> improvement, or an undesired change due to the DEM being treated
> differently between generating the RPCs and using them, I'm not sure!
>
> Daniel
>
> On Thu, 18 Apr 2024 at 15:39, Joseph McGlinchy <jmcglinchy at hydrosat.com>
> wrote:
>
> Hi Daniel,
>
> Thanks for trying with the data! I appreciate that. The gdal command I
> provided in the sample-gdal-call.txt file is the only one I could get to
> have "good" looking results. What command did you use to include the DEM?
>
> If there has been a fix between 3.6 and 3.8.3 that helps with this, then
> I'm happy to try it. I'll give that a look too.
>
> thanks!
> Joe
> ------------------------------
> *From:* Daniel Evans <daniel.fred.evans at gmail.com>
> *Sent:* Thursday, April 18, 2024 5:44 AM
> *To:* Joseph McGlinchy <jmcglinchy at hydrosat.com>
> *Cc:* Even Rouault <even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
> <gdal-dev at lists.osgeo.org>
> *Subject:* Re: [gdal-dev] gdalwarp with RPCs causing bogus output
>
> Hi Joe,
>
> Trying your commands with GDAL 3.8.3, the results don't seem to be as bad
> as what you're reporting - I don't see big strips missing. Adding in the
> DEM file as Even mentioned, it seems to be doing what I'd expect (the
> correction isn't "right" because of the georeferencing offset, but it's
> doing what it's told using the underlying DEM).
>
> Any possibility of trying a newer GDAL version in case there has been a
> fix since 3.6.0?
>
> Cheers,
> Daniel
>
> On Thu, 18 Apr 2024 at 02:35, Joseph McGlinchy via gdal-dev <
> gdal-dev at lists.osgeo.org> wrote:
>
> Hi Even,
>
> Thank you for this great response. I'll try out some of those things you
> mention on the CRS of the DEM and correcting for the geoid; I had also been
> reading that any height values used in this process need to be heights
> above the ellipsoid.
>
> I had seen good outputs from gdaltransform, both forward and inverse,
> and can get back to the corner coordinates, center coordinates, etc. in
> both lon/lat and sample/line so it could be a subtlety of the DEM + CRS, or
> the more challenging issue of stable RPCs.
>
> On the topic of stable RPCs, do you happen to know if it is enough to zero
> out the higher-order terms after the full set of coefficients are fit?
> Since there are squared terms in the higher-order coefficients, that could
> definitely be part of the issue, even if it works mathematically.
>
> thanks!
> Joe
> ------------------------------
> *From:* Even Rouault <even.rouault at spatialys.com>
> *Sent:* Wednesday, April 17, 2024 4:20 PM
> *To:* Joseph McGlinchy <jmcglinchy at hydrosat.com>; gdal-dev at lists.osgeo.org
> <gdal-dev at lists.osgeo.org>
> *Subject:* Re: [gdal-dev] gdalwarp with RPCs causing bogus output
>
>
> Hi,
>
>
> I haven't done myself that exercice but I know that computing RPC values
> that are stable enough might be challenging, so perhaps the issue is
> related to that
>
>
> A few other remarks:
>
> - your gdalwarp command line refers to  RPC_DEM_SRS and
> RPC_DEM_MISSING_VALUE but doesn't include a RPC_DEM itself, hence those are
> likely to be non effective
>
> - you generally don't want to use both RPC_HEIGHT and RPC_DEM.  RPC_HEIGHT
> is essentially useful when you don't have a DEM available, and thus
> fallback taking an average elevation
>
> - there's a subtlety regarding DEM. RPC reference for altitudes is WGS84
> ellipsoidal height, not orthometric/MSL altitude. But DEM values use
> orthometric/MSL altitude, hence a geoid correction must be applied. So
> you'd rather want to use RPC_DEM_SRS=EPSG:4326+5773 for example to use the
> EGM96 geoid (cf https://github.com/OSGeo/gdal/issues/3298). That said,
> misusing orthometric altitude vs ellipsoidal height generally accounts for
> small shifts, not big ones
>
> - you could use gdaltransform with all your -to parameter and
> input_with_RPCs.tiff to check at least that the forward RPC transformation
> path works correctly (the one from longitude, latitude -> column, line).
> And with -i to check the inverse RPC transformer. The inverse RPC
> transformer can have a hard time converging in montainous areas or for
> images off-nadir
>
>
> Even
>
>
> Le 17/04/2024 à 23:30, Joseph McGlinchy via gdal-dev a écrit :
>
> Hello,
>
> I am attempting to implement georegistration through RPC. I have the
> following information I've used to calibrate the RPC coefficients, using
> all terms for numerator and denominator for both sample and line equations.
>
>
>    - image grid, stored in tiff format with no geo-information associated
>    with it, so it reflects the imaging orientation
>    - a 'ground' grid, which corresponds to the longitude/latitude
>    coordinates for each pixel determined from a line-of-sight vector and
>    imaging system coordinates where the pseudo-rays intersect the WGS84 geoid
>    - a random sample of up to 4000 image coordinate / object space
>    coordinate pairs (I arrived at this number through trial and error; using
>    all pixels explodes RAM)
>    - elevation extracted at the longitude/latitude coordinates from a
>    DEM, in this case, SRTM, but have also tried using NASADEM
>
>
> I am able to write out an image to EPSG:4326 by populating the RPC
> metadata and using the non-georeferenced image data. However, I am
> struggling to use gdalwarp in a reliable way to orthorectify that data,
> let alone writing it out in a way that 'bakes in' the georegistration with
> the RPCs whether that is in EPSG:4326 or the local UTM zone. I see strips
> of the image "ripped out" with odd curves in various places throughout the
> image.
>
> The only way I've been able to use gdalwarp to write the image at all is
> with the following parameters (any DEM reference is to the SRTM DEM):
>
> gdalwarp --config CPL_DEBUG ON -t_srs EPSG:32610 -rpc -to
> "RPC_DEM_SRS=+proj=longlat +datum=WGS84 +no_def" -to "RPC_HEIGHT=350" -to
> "RPC_DEM_MISSING_VALUE=0.001" -to "RPC_FOOTPRINT='POLYGON ((list of polygon
> coordinates comprising the long/lat grid))'" -to "RPC_MAX_ITERATIONS=101"
> input_with_RPCs.tiff output.tiff
>
> This is the only configuration i can use to run gdalwarp successfully.
> removing any single RPC_X tranformer option gives me bogus output. The
> RPC_HEIGHT value i specify above is not close to the mean or median
> elevation of the extent of my data; mean is ~195m and median is ~150m.
>
> With the debug turned on, any other set of parameters gives me failed RPC
> convergence on several points. I am able to reproduce this regularly by
> specifying RPC_DEM=dem.tif, where dem.tif is the same data I used to
> extract elevation values when calibrating the RPCs. I am seeing normalized
> latitude and longitude values with magnitude > 1 (I checked every location
> in the image, based on the metadata, the range is not outside of [-1,1]),
> as well as normalized altitude values with magnitude > 1 (there are some,
> not many, that have magnitude of 1.75).
>
> My workflow can be summarized as:
>
>
>    1. load grids (image data, longitude, latitude)
>    2. randomly sample up to 4000 points in image coordinates, object
>    coordinates
>    1. assign z-value from SRTM DEM
>       2. evaluate if any of the points are in NODATA areas of SRTM (image
>       is coastal, so there are NODATA areas for SRTM here), if so, remove those
>       and generate more points
>       3. normalize coordinates of grids to be in [-1,1], recording
>    offsets and scale
>    4.  calibrate RPC coefficients using all terms
>    5. write out GeoTIFF with image grid for pixels, along with RPC
>    required metadata fields and CRS EPSG:4326
>
>
> System information (please let me know if more is needed)
> OS: Ubuntu 20.04 LTS (GNU/Linux 5.10.16.3-microsoft-standard-WSL2 x86_64)
> (Windows Subsystem for Linux)
> GDAL 3.6.0 (python)
>
> Thank you in advance for any insight into this process! I am happy to
> package up any of the data I am using, as well. I have placed the initial
> data from the end of Step 5 described above, along with some additional
> files, a sample gdalwarp call, and a file-list.txt, at
> https://drive.google.com/drive/folders/1BfevhKQa4ZHi_OQfiX_rUk2sqeoNVhyM?usp=sharing
>
> If more is needed for anyone interested in having a look, please let me
> know and I'll upload.
>
>
> Cheers,
> Joe
>
>
> _______________________________________________
> gdal-dev mailing listgdal-dev at lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240422/2f8db95e/attachment-0001.htm>