4

I have a raster image (photograph) that was taken from a helicopter. Along with that raster image, the system also tells me the positions of the four corners of the image in a different file. I have imported the image into R but it does not have any georeferencing builtin to it. I would like to find a way to tell R what each of the corners are so that I can plot it on a map and perform different types of spatial analysis work on the image.

Most posts on this subject specify that all I need to do is tell R about the bounding box and everything will be good to go. The issue there is that the image can have an arbitrary orientation and is not necessarily parallel to, well a parallel.

I am sure that there is an easy way to do this but I am very new to this type of work.

1

I think you would be far better off georeferencing your image in QGIS. R can do a lot of things, but its not yet a full replacement GIS. The georeferencing system in QGIS is very simple to use, free, and you can throw in some streaming satellite imagery and/or other free map data like OpenStreetMap in the background to help you position the image as correctly as possible. Your description doesn't sound like a proper airphoto, so expect to have to use a lot of control points, especially towards the edges of the image. Make sure the coordinates you were given are used for the first four control points, and then add as many as you need. Experiment with the transformation algorithms, more become available as you add more points.

1

I' think it is not a matter of R ..only. Here a small sketch howto realize the referencing stuff under Linux. You can use GDAL to create a georeferenced axis parallel image from your heli position data, assuming, that your image has the dimension width and height and the bounding box is given by may be GPS-coodinates stored in UTM-33 shown in this table:

 image-x  image-y     world-x   world-y      
 --------------------------------------
     0       0         ulx       uly
     0      height     llx       lly
   width     0         urx       ury
   width    height     lrx       lry

where variable names are:

  ulx - upper left corner x
  uly - upper left corner y
  llx - lower left corner x
  lly - lower left corner y
  urx - upper right corner x
  ury - upper right corner y
  lrx - lower right corner x
  lry - lower right corner y

of the in the world turned and shifted bounding box. With GDAL you need two steps to get the image into an axis parallel form.

  1. Create a new image with pass point and coordinate system infos using gdal_translate.

  2. Recalculate the axis parallel image with gdalwarp.

Here a small bash pseudo code:

   sid=32633 # Datum for UTM-33

   hg=<IMAGE_HEIGHT>
   wd=<IMAGE_WIDTH>
   ulx=<corner-x from gps>
   uly=<corner-y from gps>       
   ...

   gsd=10 # Target resolution lets say 10m

   # 1. Step
   gdal_translate -of GTiff\
    -gcp   0     0 $ulx   $uly \
    -gcp   0   $hg $llx   $lly \
    -gcp $wd     0 $urx   $ury \
    -gcp $wd   $hg $lrx   $lry \
    -a_srs epsg:$sid\
     in.tif out-pp.tif

    # 2. Step resampling
    gdalwarp -of GTiff  -r bilinear \
             -tr $gsd out-pp.tif out-ref.tif

With the tool gdalinfo you can read the info tags of the new georeferenced axis parallel image. An example

  $ gdalinfo out-ref.tif 
  Driver: GTiff/GeoTIFF
  Files: out-ref.tif
  Size is 5001, 6868
  Coordinate System is:
  PROJCS["WGS 84 / UTM zone 33N",
      ...
     AUTHORITY["EPSG","32633"]]
  Origin = (370000.000000000000000,6068670.000000000000000)
  Pixel Size = (9.998000399920016,-9.998543972044263)
  ...

From here you get the info offsX/Y and pixSizeX/Y, as you see it is not exaclty 10m scaled. You can switch between the both worlds (axis parallel image and real world) using the following functions and go on in R.

  calcWorldPos <- function (pixX ,pixY, offsX, offsY, pixSizeX, pixSizeY)

     X = offsX + pixX * pixSizeX;
     Y = offsY - pixY * pixSizeY;

     return c(x,y);
  }

  calcPixPos <- function( worldX, worldY, offsX, offsY, pixSizeX, pixSizeY) {

     pixX = floor(( worldX - offsX ) / pixSizeX );
     pixY = floor(( offsY - worldY ) / pixSizeY );

     return c(pixX, pixY); 
  }

There is also a package rgdal to work with georeferenced images in a more closed form and without leaving R but I miss the convenient tools gdalinfo, gdal_translate and gdalwarp.

  • For some images, this approach might sort of work. It cannot succeed all the time, though, because it attempts to approximate a projective transformation (which requires eight real parameters) by an affine one (which uses only six parameters). In particular, mages with wide fields of view require substantial internal "undistortion" that won't be accomplished by this approach. – whuber Apr 10 '15 at 18:09
  • Can you help me understand what a wide FOV might be? The altitude this image was taken at was fairly low. – Justace Clutter Apr 10 '15 at 18:29
  • Huck, Thanks for the informative answer. As I am really new to this I will need to take some time to digest this. I do not know what the physical width and height are of my images. All I have is an image that is 2048x2048 and then four lat/lon points for the image edge coordinates. It seems as though I need to know what the image width is... Wait, can't I get that from just calculating the distance between the four given lat/lon coordinates? Also can you give some guidance on how to select the appropriate resolution? Thanks for your inputs again. – Justace Clutter Apr 11 '15 at 0:25
  • @whuber Yes your comment is quite correct, but to get the distortion stuff inside the rectification of the image, we need the camera parameters, position (may be of the corners) and the IMU data of the platform (omega, kappa and phi) and the tilt of the surface..., and..., which are not present. For a first step into the stuff, I think it's worth to choose a regional common sense used projection, put the gps lon/lat data into it and work on.. – huckfinn Apr 11 '15 at 18:46
  • THe resolution of the image is given by the cell size of the sensor, the focal length and the flight height of your platform (and the camera constant) ..see en.wikipedia.org/wiki/Lens_%28optics%29#Imaging_properties. – huckfinn Apr 11 '15 at 18:55

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.