3

Now I have pleiades raster (img.tif) and a RPC model (rpc.XML). In the rpc.XML, there are the following lines which can be used to reset geo information.

  <RFM_Validity>
    <Direct_Model_Validity_Domain>
      <FIRST_ROW>1</FIRST_ROW>
      <FIRST_COL>1</FIRST_COL>
      <LAST_ROW>19676</LAST_ROW>
      <LAST_COL>18665</LAST_COL>
    </Direct_Model_Validity_Domain>
    <Inverse_Model_Validity_Domain>
      <FIRST_LON>7.393918390116374</FIRST_LON>
      <FIRST_LAT>46.89414037853461</FIRST_LAT>
      <LAST_LON>7.531577540047203</LAST_LON>
      <LAST_LAT>46.98982504658424</LAST_LAT>
    </Inverse_Model_Validity_Domain>
    <LONG_SCALE>0.06584113709872685</LONG_SCALE>
    <LONG_OFF>7.462641690331035</LONG_OFF>
    <LAT_SCALE>0.04690424904394419</LAT_SCALE>
    <LAT_OFF>46.94198271255942</LAT_OFF>
    <HEIGHT_SCALE>245</HEIGHT_SCALE>
    <HEIGHT_OFF>765</HEIGHT_OFF>
    <SAMP_SCALE>9332</SAMP_SCALE>
    <SAMP_OFF>9333</SAMP_OFF>
    <LINE_SCALE>9837.499999999996</LINE_SCALE>
    <LINE_OFF>9838.499999999996</LINE_OFF>
  </RFM_Validity>

And I use the following code to reset geo info of img.tif with the value below.
xmin = 7.393918390116374
ymin = 46.89414037853461
xmax = 7.531577540047203
ymax = 46.98982504658424
first_col = 1
last_col = 18665
first_row = 1
last_row = 19676

def set_geo_info(in_path, out_path, xmin, ymin, xmax, ymax,
                 first_row, first_col, last_row, last_col, *args):
  '''
  xmin, ymin, xmax, ymax are the longitude and latitude.
  '''
  print('Set Geo Information ...')
  xres = (xmax - xmin) / (last_col - first_col + 1)
  yres = (ymax - ymin) / (last_row - first_row + 1)

  geotransform = (xmin, xres, 0, ymax, 0, -yres)

  ds = gdal.Open(in_path)
  band = ds.GetRasterBand(1)
  arr = band.ReadAsArray()
  [y, x] = arr.shape

  driver = gdal.GetDriverByName("GTiff")
  outdata = driver.Create(out_path, x, y, 1, gdal.GDT_UInt16)
  outdata.SetGeoTransform(geotransform)
  outdata.GetRasterBand(1).WriteArray(arr)
  outdata.FlushCache()

  outdata = None
  arr = None
  band = None
  ds = None

  print("Done.")

The result is kind of correct, but not entirely.

Position shifting wrt shapefile (WGS84)

As you can see, the red line is a shapefile (EPSG21781), while the railway in the result is not entirely overlap with shapefile. The distance of shifting is not consistent for this whole image. The displacement is larger near the edge of this image.

My questions:

  1. Am I heading the right direction to reset geo information?
  2. If 1. is correct, I guess the value which I fed into the gdal is wrong. How do I find the correct value from RPC model (rpc.XML)?
  3. If 1. and 2. is wrong, how do I do it in a right way?
  • Why downvote? Please criticize so I can improve my asking skill. – Lion Lai Jan 18 '19 at 18:14
1

You can warp the image using the RPCs directly with GDAL, e.g.

gdalwarp -r cubic -to RPC_DEM=my_dem.tif -rpc -t_srs epsg:4326  my_DIMAP_file.XML my_dimap_file_warped.tif

In order to get a good result it is advisable to use a high resolution DEM. You can do the same programmatically if you prefer by using gdal.Warp()

|improve this answer|||||
  • Thank you. It works, although there's still small displacement, but it's better than before. I use python elevation packge to download DEM. I think it's the reason why there is still small displacement in the output. – Lion Lai Jan 29 '19 at 13:53
  • I don't know the details of the product that you have, but it is quite possible that it is a so-called "Ortho-ready" product ( see orfeo-toolbox.org/CookBook/recipes/… ) which means the RPCs were created by projecting to a constant height instead of a DEM. In that case the positional accuracy that can be achieved would indeed be limited. – Markus M. Jan 30 '19 at 6:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.