How to use the value in rpc.XML to reproject raster.tif?

lionlai1989 commented 5 years ago

Guys: Now I have pleiades raster (raster.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 raster.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.")

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.

However, the output DSM of s2p is really consistent with railway's position. grafik As you can see, the blue railway is in between the building perfectly.

So, my question is the following.

Am I heading the right direction to reproject raster.tif?
If 1. is correct, I guess the value which I fed into the gdal is wrong. How to find the correct value from RPC model (rpc.XML).
How does s2p use rpc.XML to figure out the coordinate system?

Thank you. Any help will be appreciated.

gfacciol commented 5 years ago

Dear @lionlai1989, What you are observing is the difference between an orthophoto and an image as it comes from the satellite. Although the corners of the raster images are geolocated, this doesn't mean that the whole image is an orthographic projection. In this case, it is clear from the photo that it is not since you see the sides of the buildings. However, when you build the DSM this, model is projected on a geographic grid, so it maches your railway path. Roughly, what you need to do is to project the image on the DSM to obtain on orthoimage. I hope the explication is clear enough.

So the answers to your questions are: 1) yes, 2) nothing is wrong with how you're using the RPC, 3) yes.

Best, gabriele

lionlai1989 commented 5 years ago

@gfacciol Thank you for your explanation. Reprojection a raster doesn't mean orthorectifying a raster.

MISS3D / s2p

How to use the value in rpc.XML to reproject raster.tif? #187