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, <jmcglin...@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.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> > 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> > *Sent:* Thursday, April 18, 2024 5:44 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, > > 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> 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> > *Sent:* Wednesday, April 17, 2024 4:20 PM > *To:* Joseph McGlinchy <jmcglin...@hydrosat.com>; gdal-dev@lists.osgeo.org > <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 > 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@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@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
