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.ev...@gmail.com>
Sent: Thursday, April 18, 2024 9:03 AM
To: Joseph McGlinchy <jmcglin...@hydrosat.com>
Cc: Even Rouault <even.roua...@spatialys.com>; gdal-dev@lists.osgeo.org
<gdal-dev@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
<jmcglin...@hydrosat.com<mailto:jmcglin...@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.ev...@gmail.com<mailto:daniel.fred.ev...@gmail.com>>
Sent: Thursday, April 18, 2024 5:44 AM
To: Joseph McGlinchy <jmcglin...@hydrosat.com<mailto:jmcglin...@hydrosat.com>>
Cc: Even Rouault
<even.roua...@spatialys.com<mailto:even.roua...@spatialys.com>>;
gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
<gdal-dev@lists.osgeo.org<mailto:gdal-dev@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@lists.osgeo.org<mailto:gdal-dev@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.roua...@spatialys.com<mailto:even.roua...@spatialys.com>>
Sent: Wednesday, April 17, 2024 4:20 PM
To: Joseph McGlinchy <jmcglin...@hydrosat.com<mailto:jmcglin...@hydrosat.com>>;
gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
<gdal-dev@lists.osgeo.org<mailto:gdal-dev@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
*
assign z-value from SRTM DEM
*
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 list
gdal-dev@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
https://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@lists.osgeo.org<mailto:gdal-dev@lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
