Orthorectification digital camera

From GRASS-Wiki
Jump to: navigation, search

Geocoding of photographs taken with digital handheld camera

To geocode oblique photographs taken with a digital handheld camera, we need a digital elevation model (DEM) and a topographic reference map. The elevation model is required to normalize the terrain undulation. The raster resolution of the DEM should be similar to the resolution of the aerial photo. Usually such high-resolution DEMs (raster cell length below 1m) are not easily available. Therefore, interpolation of a given DEM to a higher resolution is recommended to minimize displacement effects. The topographic reference map is needed to find corresponding ground control points (GCPs) between aerial photo and this map for geocoding. The map scale of the reference map should at least match the average scale of the aerial photo (for example 1:5000), but should be probably larger to provide more details for improved geocoding accuracy. The ground control points are required to register the image xy coordinate system to a georeferenced coordinate system (like UTM or Gauss-Krüger).

Metadata like the time stamp or the focal length are taken from EXIF data (e.g. with "jhead", http://www.sentex.net/~mwandel/jhead/ ).

PROCEDURE

Ortho-rectification of oblique photographs

In projected location:


- get approximate camera position:

  d.where/d.what.rast:
              EAST:             NORTH:
   664194.71461534   5104463.03343466

 height: 450

  Pitch: 20
  Roll: 0
  Yaw (against north): -45

See also: http://en.wikipedia.org/wiki/Pitch_%28flight%29

######################################

g.region rast=cf83 -p
projection: 0 (x,y)
zone:       0
north:      798
south:      0
west:       0
east:       1200
nsres:      1
ewres:      1
rows:       798
cols:       1200

 -> img size in Pixel (NOTE: newer cameras provide this in exif metadata!)

i.ortho.photo

Menu 5. Compute image-to-photo transformation

Some considerations:

 Calibrated Focal Length mm.: guessing 50 for a digital handheld camera (airborne cameras have often a focal length of 150 mm)

 Fiducial marks (positive values, use g.region -p):

                        y
                        |        2
                        +--------+--------+
                        |                 |
                        |                 |
                        |                 |
                      1 +                 + 3
                        |                 |
                        |                 |
                        |                 |
                        +--------+--------+------> x
                      (0,0)      4

  also: xmax = 1200/2 = 600
        ymax = 798/2  = 399

Digital cameras have a resolution
Fiducial marks have to be given in millimeters (so, to calculate)
         and not in pixel:

first calc. DPI, then derive "pseudo''millimeters from that
for Fiducial marks:

E.g.: Lenght of an object (eg bridge) in pixel (use d.what.rast + Pytaghoras)
compared to length of bridge in meters on reference map (d.measure).

  eg 226 pixel versus 75 meter in UTM map

 [ DPI = pixel/inch = pixel / 2.54cm ]

->  226/7500 * 2.54cm = 0.07653867 DPI (check)

Since we calculate in *millimeter* and not DPI, we calc. for Fiducial marks:
  226/7500 * 0.1 = 0.003013333

Multiply all pixel with this value:

 1200 * 0.003013333 = 3.616000
 1200 * 0.003013333 /2 = 1.808000
 798 * 0.003013333 = 2.404640
 798 * 0.003013333 /2 = 1.202320

Ergo:
             Fiducial marks:
             1       1_____      0_________      1.2______
             2       2_____      1.8______       2.4______
             3       3_____      3.61_____       1.2______
             4       4_____      3.61______      0_________

(be a bit more precise...)

##################
Then use menu 5 in i.ortho.photo. It's a challenge to hit the center
of the image borders... try smallest RMS


Menu 6. Initialize exposure station parameters

     Initial Camera Exposure X-coordinate Meters:           664194________
     Initial Camera Exposure Y-coordinate Meters:           5104463_______
     Initial Camera Exposure Z-coordinate Meters:           200___________
     Initial Camera Omega (roll) degrees:                   0_____________
     Initial Camera Phi  (pitch) degrees:                   20____________
     Initial Camera Kappa  (yaw) degrees:                   -45___________

Menu 7. Compute ortho-rectification parameters

Use mouse now

Search ground control points, minimize RMS (see Analyse menu entry) to 
ideally half a pixel.

Menu 8. Ortho-rectify imagery files

Enter an extension to be appended to rectified maps:
.ortho______________


Compute local camera angle? (y/n) [n] y

Enter a name for the camera angle map:
22_1684.camera_angle__________


Overwrite maps in target location/mapset? (y/n) [n] y


Please select one of the following options
 1. Use the current window in the target location
 2. Determine the smallest window which covers the image
> 2

Enter desired target resolution, or
 RETURN   to determine it automatically:
> 2.0

Please select one of the following interpolation methods
 1. nearest neighbor
 2. bilinear
 3. bicubic
 4. bilinear with fallback
 5. bicubic with fallback
> 2


Enter amount of memory to use in MB, or
RETURN   use 100 MB
> 1000

Now the orthophoto is generated and will be stored in the target location.

Hints: How to calculate Fiducial marks, if no DPI known


- Pixel extent: g.region rast=oblique

                        y
                        |        2
                        +--------+--------+
                        |                 |
                        |                 |
                        |                 |
                      1 +                 + 3
                        |                 |
                        |                 |
                        |                 | 
                        +--------+--------+------> x
                      (0,0)      4

e.g.
        xmax = 1200/2 = 600
        ymax = 798/2  = 399

Important: while digital handheld cameras provide DPI and chip size, one has to estimate everything if it is unknown. Fiducial marks have to be given in millimeters, not in pixel!

Calculation of Pseudomillimeter:

- measure in XY location length of an object in pixel (with d.what.rast + 
   Pytaghoras)
- measure in UTM location length of same object in reference map (d.measure)
   in meters

Example: bridge 226 pixel  versus 75 meter in UTM.

[ DPI = pixel/inch = pixel / 2.54cm ]

->  226/7500 * 2.54cm = 0.07653867 DPI

BUT: We want *millimeter* and not DPI for the fiducial marks:
226/7500 * 0.1 = 0.003013333

Finally multiply pixels with this factor:

1200 * 0.003013333 = 3.616000
1200 * 0.003013333 /2 = 1.808000
798 * 0.003013333 = 2.404640
798 * 0.003013333 /2 = 1.202320

Ergo:

             Fiducial marks:
             1       1_____      0_________      1.2______
             2       2_____      1.8______       2.4______
             3       3_____      3.61_____       1.2______
             4       4_____      3.61______      0_________


Suggestions: 
 - take long objects,
 - average over 2 longitudinal and transversal oriented objects

See also