Orthorectification WorldView-3 Error

gis123 · September 11, 2020, 7:28pm

Hi All,

I’m trying to complete orthorectification on my WV-3 satellite imagery, but after I run it on my ortho-ready imagery, literally nothing in the image changes. I’m comparing the “orthorectified” product to basemaps and OSM and the images are probably about 20m off. I’ve played around with all the OTB parameters and nothing seems to work.

I have checked and I do have my RPCs stored in an .RPB file and I am using the SRTM DEM which is about ~30m in resolution. I’ve verified that the DEM is in the same projection system as the satellite image I am trying to orthorectify. I’ve also tried using the GEOID file that is available for download and that doesn’t work either.

I’ve tried attempting this in QGIS (3.8.1) with OTB and also with MonteVerdi (version OTB-6.6.1-Win64), to no avail. I am using Windows 10. Right now I’m at a loss how to correct this data and haven’t been able to find other answers to similar issues using this platform.

I appreciate any response!

yannick · September 14, 2020, 3:18pm

Hi,

Could you describe more precisely how you used OTB (application and parameters used) ?
I suppose you have used the OrthoRectification application : the SRTM DEM and the EGM 96 geoid are fine to process ortho-rectification.
To be sure that OTB can open and process your WV-3 product, you could use either the ReadImageInfo application, or gdalinfo :

otbcli_ReadImageInfo -in yourImage.tif
gdalinfo yourImage.tif

Be careful that in some viewers (QGIS, monteverdi), the images are displayed with a default projection system : the offset you may observe may be considered with respect to that default projection.

gdalinfo will display corner coordinates only if the image had been georeferenced : can you confirm that your resulting image had no geographic coordinates ?

Hope that helps

Yannick

gis123 · September 15, 2020, 5:00pm

Hi Yannick,

Thanks for the reply. The main issue is that the software returns no errors and finishes the orthorectification but the pre-ortho and post-ortho image look the same. It doesn’t seem to be properly utilizing the DEM. When I use gdalinfo on the tiff before ortho, it results in this message:

Size is 27902, 22369
Coordinate System is:
PROJCRS[“WGS 84 / UTM zone 16N”,
BASEGEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223563,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]],
CONVERSION[“UTM zone 16N”,
METHOD[“Transverse Mercator”,
ID[“EPSG”,9807]],
PARAMETER[“Latitude of natural origin”,0,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8801]],
PARAMETER[“Longitude of natural origin”,-87,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8802]],
PARAMETER[“Scale factor at natural origin”,0.9996,
SCALEUNIT[“unity”,1],
ID[“EPSG”,8805]],
PARAMETER[“False easting”,500000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8806]],
PARAMETER[“False northing”,0,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT[“metre”,1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT[“metre”,1]],
USAGE[
SCOPE[“unknown”],
AREA[“World - N hemisphere - 90┬░W to 84┬░W - by country”],
BBOX[0,-90,84,-84]],
ID[“EPSG”,32616]]
Data axis to CRS axis mapping: 1,2
Origin = (666468.900000000023283,3881626.799999999813735)
Pixel Size = (0.300000000000000,-0.300000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 666468.900, 3881626.800) ( 85d10’27.87"W, 35d 3’49.18"N)
Lower Left ( 666468.900, 3874916.100) ( 85d10’32.71"W, 35d 0’11.45"N)
Upper Right ( 674839.500, 3881626.800) ( 85d 4’57.56"W, 35d 3’44.08"N)
Lower Right ( 674839.500, 3874916.100) ( 85d 5’ 2.64"W, 35d 0’ 6.36"N)
Center ( 670654.200, 3878271.450) ( 85d 7’45.19"W, 35d 1’57.80"N)
Band 1 Block=27902x1 Type=UInt16, ColorInterp=Gray
Overviews: 13951x11185, 6976x5593
Band 2 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593
Band 3 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593
Band 4 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593

After using the orthorectification OTB tool, gdalinfo reads this:
Size is 27902, 22369
Coordinate System is:
PROJCRS[“WGS 84 / UTM zone 16N”,
BASEGEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223563,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]],
CONVERSION[“UTM zone 16N”,
METHOD[“Transverse Mercator”,
ID[“EPSG”,9807]],
PARAMETER[“Latitude of natural origin”,0,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8801]],
PARAMETER[“Longitude of natural origin”,-87,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8802]],
PARAMETER[“Scale factor at natural origin”,0.9996,
SCALEUNIT[“unity”,1],
ID[“EPSG”,8805]],
PARAMETER[“False easting”,500000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8806]],
PARAMETER[“False northing”,0,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT[“metre”,1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT[“metre”,1]],
USAGE[
SCOPE[“unknown”],
AREA[“World - N hemisphere - 90┬░W to 84┬░W - by country”],
BBOX[0,-90,84,-84]],
ID[“EPSG”,32616]]
Data axis to CRS axis mapping: 1,2
Origin = (666468.900000000023283,3881626.799999999813735)
Pixel Size = (0.300000000000000,-0.300000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 666468.900, 3881626.800) ( 85d10’27.87"W, 35d 3’49.18"N)
Lower Left ( 666468.900, 3874916.100) ( 85d10’32.71"W, 35d 0’11.45"N)
Upper Right ( 674839.500, 3881626.800) ( 85d 4’57.56"W, 35d 3’44.08"N)
Lower Right ( 674839.500, 3874916.100) ( 85d 5’ 2.64"W, 35d 0’ 6.36"N)
Center ( 670654.200, 3878271.450) ( 85d 7’45.19"W, 35d 1’57.80"N)
Band 1 Block=27902x1 Type=UInt16, ColorInterp=Gray
Overviews: 13951x11185, 6976x5593
Band 2 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593
Band 3 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593
Band 4 Block=27902x1 Type=UInt16, ColorInterp=Undefined
Overviews: 13951x11185, 6976x5593

(spark) D:\imagery\apollo-chattanooga-wv3>gdalinfo 20APR18163440-S2AS_O1-013028683010_01_P001.tif
Driver: GTiff/GeoTIFF
Files: 20APR18163440-S2AS_O1-013028683010_01_P001.tif
Size is 27984, 22466
Coordinate System is:
PROJCRS[“WGS 84 / UTM zone 16N”,
BASEGEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223563,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]],
CONVERSION[“UTM zone 16N”,
METHOD[“Transverse Mercator”,
ID[“EPSG”,9807]],
PARAMETER[“Latitude of natural origin”,0,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8801]],
PARAMETER[“Longitude of natural origin”,-87,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8802]],
PARAMETER[“Scale factor at natural origin”,0.9996,
SCALEUNIT[“unity”,1],
ID[“EPSG”,8805]],
PARAMETER[“False easting”,500000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8806]],
PARAMETER[“False northing”,0,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT[“metre”,1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT[“metre”,1]],
USAGE[
SCOPE[“unknown”],
AREA[“World - N hemisphere - 90┬░W to 84┬░W - by country”],
BBOX[0,-90,84,-84]],
ID[“EPSG”,32616]]
Data axis to CRS axis mapping: 1,2
Origin = (666456.900000000023283,3881641.799999999813735)
Pixel Size = (0.299999999999998,-0.299999999999992)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 666456.900, 3881641.800) ( 85d10’28.33"W, 35d 3’49.68"N)
Lower Left ( 666456.900, 3874902.000) ( 85d10’33.19"W, 35d 0’11.00"N)
Upper Right ( 674852.100, 3881641.800) ( 85d 4’57.05"W, 35d 3’44.56"N)
Lower Right ( 674852.100, 3874902.000) ( 85d 5’ 2.16"W, 35d 0’ 5.89"N)
Center ( 670654.500, 3878271.900) ( 85d 7’45.18"W, 35d 1’57.81"N)
Band 1 Block=27984x1 Type=UInt16, ColorInterp=Gray
NoData Value=0
Band 2 Block=27984x1 Type=UInt16, ColorInterp=Undefined
NoData Value=0
Band 3 Block=27984x1 Type=UInt16, ColorInterp=Undefined
NoData Value=0
Band 4 Block=27984x1 Type=UInt16, ColorInterp=Undefined
NoData Value=0

yannick · September 16, 2020, 12:21pm

Hello,

According to me, your input image is already orthorectified : its projection system is
WGS 84 / UTM 16N, and gdalinfo displays the corners coordinates (lat/lon).
In an image with a raw geometry, the corner coordinates would be only a pixel index.

Your output image seems to have exactly the same characteristics : you can use the CompareImages to check if they are the same :

otbcli_CompareImages -ref.in -meas.in
It will display some stats and indicate if they are identical.

Maybe we should improve OrthoRectification application to prevent being used on an already georeferenced image.

Yannick