Main Content

Georeferencing an Image to an Orthotile Base Layer

This example shows how to register an image to an earth coordinate system and create a new "georeferenced" image. It requires Image Processing Toolbox™ in addition to Mapping Toolbox™.

In this example, all georeferenced data are in the same earth coordinate system, the Massachusetts State Plane system (using the North American Datum of 1983 in units of meters). This defines our "map coordinates." The georeferenced data include an orthoimage base layer and a vector road layer.

The data set to be georeferenced is a digital aerial photograph covering parts of the village of West Concord, Massachusetts, collected in early spring, 1997.

Step 1: Render Orthoimage Base Tiles in Map Coordinates

The orthoimage base layer is structured into 4000-by-4000 pixel tiles, with each pixel covering exactly one square meter in map coordinates. Each tile is stored as a TIFF image with a world file. The aerial photograph of West Concord overlaps two tiles in the orthoimage base layer. (For convenience, this example actually works with two 2000-by-2000 sub-tiles extracted from the larger 4000-by-4000 originals. They have the same pixel size, but cover only the area of interest.)

Read the two orthophoto base tiles required to cover the extent of the aerial image.

[baseImage1,R1,cmap1] = readgeoraster("concord_ortho_w.tif");
[baseImage2,R2,cmap2] = readgeoraster("concord_ortho_e.tif");

Display the images in their correct spatial positions.

mapshow(baseImage1,cmap1,R1)
ax1 = gca;
mapshow(ax1,baseImage2,cmap2,R2)
title("Map View, Massachusetts State Plane Coordinates");
xlabel("Easting (meters)");
ylabel("Northing (meters)");

Step 2: Register Aerial Photograph to Map Coordinates

This step uses functions cpselect, cpstruct2pairs, fitgeotform2d, and imwarp, and method projtform2d/transformPointsForward, from the Image Processing Toolbox together with map raster reference objects from Mapping Toolbox. Together, they enable georegistration based on control point pairs that relate the aerial photograph to the orthoimage base layer.

Read the aerial photo.

inputImage = imread("concord_aerial_sw.jpg");
figure
imshow(inputImage)
title("Unregistered Aerial Photograph")

Both orthophoto images are indexed images but cpselect only takes grayscale images, so convert them to grayscale.

baseGray1 = im2uint8(ind2gray(baseImage1,cmap1));
baseGray2 = im2uint8(ind2gray(baseImage2,cmap2));

Downsample the images by a factor of 2, then pick two separate sets of control point pairs: one for points in the aerial image that appear in the first tile, and another for points that appear in the second tile. Save the control point pairs to the base workspace as control point structures named cpstruct1 and cpstruct2.

n = 2; % downsample by a factor n
load mapexreg.mat % load some points that were already picked

Optionally, edit or add to the pre-picked points using cpselect. Note that you need to work separately on the control points for each orthotile.

cpselect(inputImage(1:n:end,1:n:end,1),...
         baseGray1(1:n:end,1:n:end),cpstruct1);
cpselect(inputImage(1:n:end,1:n:end,1),...
         baseGray2(1:n:end,1:n:end),cpstruct2);

This tool helps you pick pairs of corresponding control points. Control points are landmarks that you can find in both images, like a road intersection, or a natural feature. Three pairs of control points have already been picked for each orthophoto tile. If you want to proceed with these points, go to Step 3: Infer and apply geometric transformation. If you want to add some additional pairs of points, do so by identifying landmarks and clicking on the images. Save control points by choosing the File menu, then the Save Points to Workspace option. Save the points, overwriting variables cpstruct1 and cpstruct2.

Step 3: Infer and Apply Geometric Transformation

Extract control point pairs from the control point structures.

[input1,base1] = cpstruct2pairs(cpstruct1);
[input2,base2] = cpstruct2pairs(cpstruct2);

Account for downsampling by factor n.

input1 = n*input1 - 1;
input2 = n*input2 - 1;
base1 = n*base1 - 1;
base2 = n*base2 - 1;

Transform base image coordinates into map (State Plane) coordinates.

[x1, y1] = intrinsicToWorld(R1,base1(:,1),base1(:,2));
[x2, y2] = intrinsicToWorld(R2,base2(:,1),base2(:,2));

Combine the two sets of control points and infer a projective transformation. (The projective transformation should be a reasonable choice, since the aerial image is from a frame camera and the terrain in this area is fairly gentle.)

input = [input1; input2];
spatial = [x1 y1; x2 y2];

tform = fitgeotform2d(input,spatial,"projective")
tform = 

  projtform2d with properties:

    Dimensionality: 2
                 A: [3×3 double]

Compute and plot (2D) residuals.

residuals = transformPointsForward(tform, input) - spatial;
figure
plot(residuals(:,1),residuals(:,2),"+")
title("Control Point Residuals");
xlabel("Easting offset (meters)");
ylabel("Northing offset (meters)");
xlim([-4 4])
ylim([-4 4])
axis equal

Predict corner locations for the registered image, in map coordinates, and connect them to show the image outline.

mInput = size(inputImage,1);
nInput = size(inputImage,2);

inputCorners = 0.5 ...
    + [0        0;
       0        mInput;
       nInput   mInput;
       nInput   0;
       0        0];

outputCornersSpatial = transformPointsForward(tform, inputCorners);

outputCornersX = outputCornersSpatial(:,1);
outputCornersY = outputCornersSpatial(:,2);

figure(ax1.Parent)
line(outputCornersX,outputCornersY,"Color","w")

Calculate the average pixel size of the input image (in map units).

pixelSize = [hypot( ...
    outputCornersX(2) - outputCornersX(1), ...
    outputCornersY(2) - outputCornersY(1)) / mInput, ...
 hypot( ...
    outputCornersX(4) - outputCornersX(5), ...
    outputCornersY(4) - outputCornersY(5)) / nInput]
pixelSize =

    0.9096    0.8905

Variable pixelSize gives a starting point from which to select a width and height for the pixels in our georegistered output image. In principle we could select any size at all for our output pixels. However, if we make them too small we waste memory and computation time, ending up with a big (many rows and columns) blurry image. If we make them too big, we risk aliasing the image as well as needlessly discarding information in the original image. In general we also want our pixels to be square. A reasonable heuristic is to select a pixel size that is slightly larger than max(pixelSize) and is a "reasonable" number (i.e., 0.95 or 1.0 rather than 0.9096). Here we chose 1, which means that each pixel in our georegistered image will cover one square meter on the ground. It's "nice" to have georegistered images that align along integer map coordinates for display and analysis.

outputPixelSize = 1;

Choose world limits that are integer multiples of the output pixel size.

xWorldLimits = outputPixelSize ...
    * [floor(min(outputCornersX) / outputPixelSize), ...
        ceil(max(outputCornersX) / outputPixelSize)]

yWorldLimits = outputPixelSize ...
    * [floor(min(outputCornersY) / outputPixelSize), ...
        ceil(max(outputCornersY) / outputPixelSize)]
xWorldLimits =

      208154      209693


yWorldLimits =

      911435      912583

Display a bounding box for the registered image.

line(xWorldLimits([1 1 2 2 1]),yWorldLimits([2 1 1 2 2]),"Color","red")

The dimensions of the registered image will be as follows:

mOutput = diff(yWorldLimits) / outputPixelSize
nOutput = diff(xWorldLimits) / outputPixelSize
mOutput =

        1148


nOutput =

        1539

Create an Image Processing Toolbox referencing object for the registered image.

R = imref2d([mOutput nOutput],xWorldLimits,yWorldLimits)
R = 

  imref2d with properties:

           XWorldLimits: [208154 209693]
           YWorldLimits: [911435 912583]
              ImageSize: [1148 1539]
    PixelExtentInWorldX: 1
    PixelExtentInWorldY: 1
    ImageExtentInWorldX: 1539
    ImageExtentInWorldY: 1148
       XIntrinsicLimits: [0.5000 1.5395e+03]
       YIntrinsicLimits: [0.5000 1.1485e+03]

Create a map raster reference object, which is the Mapping Toolbox counterpart to an Image Processing Toolbox referencing object.

Rmap = maprefcells(R.XWorldLimits,R.YWorldLimits, R.ImageSize, ...
    "ColumnsStartFrom","north")
Rmap = 

  MapCellsReference with properties:

            XWorldLimits: [208154 209693]
            YWorldLimits: [911435 912583]
              RasterSize: [1148 1539]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 1
      CellExtentInWorldY: 1
    RasterExtentInWorldX: 1539
    RasterExtentInWorldY: 1148
        XIntrinsicLimits: [0.5 1539.5]
        YIntrinsicLimits: [0.5 1148.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []


Apply the geometric transformation using imwarp. Flip its output so that the columns run from north to south.

registered = flipud(imwarp(inputImage, tform,"OutputView",R));
figure
imshow(registered)

You can write the registered image as a TIFF with a world file, thereby georeferencing it to our map coordinates.

imwrite(registered,"concord_aerial_sw_reg.tif");
worldfilewrite(Rmap,getworldfilename("concord_aerial_sw_reg.tif"));

Step 4: Display the Registered Image in Map Coordinates

Display the registered image on the same (map coordinate) axes as the orthoimage base tiles. The registered image does not completely fill its bounding box, so it includes null-filled triangles. Create an alpha data mask to make these fill areas render as transparent.

alphaData = registered(:,:,1);
alphaData(alphaData~=0) = 255;

figure
mapshow(baseImage1,cmap1,R1)
ax2 = gca;
mapshow(ax2,baseImage2,cmap2,R2)
title("Map View, Massachusetts State Plane Coordinates");
xlabel("Easting (meters)");
ylabel("Northing (meters)");

mInput = mapshow(ax2,registered,Rmap);
mInput.AlphaData = alphaData;

You can assess the registration by looking at features that span both the registered image and the orthophoto images.

Step 5: Overlay Vector Road Layer (from Shapefile)

Use shapeinfo and readgeotable to learn about the attributes of the vector road layer. Render the roads on the same axes and the base tiles and registered image.

roadsfile = "concord_roads.shp";
roadinfo = shapeinfo(roadsfile);
roads = readgeotable(roadsfile);

Create symbolization based on the CLASS attribute (the type of road). Note that very minor roads have CLASS values equal to 6, so don't display them.

roadspec = makesymbolspec("Line",{'CLASS',6,'Visible','off'});

mapshow(ax2,roads,"SymbolSpec",roadspec,"Color","cyan")

Observe that the roads line up quite well with the roads in the images. Two obvious linear features in the images are not roads but railroads. The linear feature that trends roughly east-west and crosses both base tiles is the Fitchburg Commuter Rail Line of the Massachusetts Bay Transportation Agency. The linear feature that trends roughly northwest-southeast is the former Framingham-Lowell secondary line.

Credits

concord_orthow_w.tif, concord_ortho_e.tif, concord_roads.shp:

  Office of Geographic and Environmental Information (MassGIS),
  Commonwealth of Massachusetts  Executive Office of Environmental Affairs
  http://www.state.ma.us/mgis
  For more information, run:
  >> type concord_ortho.txt
  >> type concord_roads.txt

concord_aerial_sw.jpg

  Visible color aerial photograph courtesy of mPower3/Emerge.
  For more information, run:
  >> type concord_aerial_sw.txt

See Also

Functions

Objects