Creating a Half-Resolution Georeferenced Image

Open Live Script

This example shows how to create a half-resolution version of a georeferenced TIFF image, using referencing objects and Image Processing Toolbox™ functions ind2gray and imresize.

Step 1: Import a Georeferenced TIFF Image

Read an indexed-color TIFF image and convert it to grayscale. The size of the image is 2000-by-2000.

[X, cmap] = imread('concord_ortho_w.tif');
I_orig = ind2gray(X, cmap);

Read the corresponding world file. Each image pixel covers a one-meter square on the map. 

R_orig = worldfileread('concord_ortho_w.tfw','planar',size(X));

Choose a convenient format for displaying the result.

currentFormat = get(0,'format');
format short g
R_orig
R_orig = 
  MapCellsReference with properties:

            XWorldLimits: [207000 209000]
            YWorldLimits: [911000 913000]
              RasterSize: [2000 2000]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 1
      CellExtentInWorldY: 1
    RasterExtentInWorldX: 2000
    RasterExtentInWorldY: 2000
        XIntrinsicLimits: [0.5 2000.5]
        YIntrinsicLimits: [0.5 2000.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []

Step 2: Resize the Image to Half Its Original Size

Halve the resolution, creating a smaller (1000-by-1000) image.

I_half = imresize(I_orig, size(I_orig)/2, 'bicubic');

Step 3: Construct a Referencing Object for the Resized Image

The resized image has the same limits as the original, just a smaller size, so copy the referencing object and reset its RasterSize property.

R_half = R_orig;
R_half.RasterSize = size(I_half)
R_half = 
  MapCellsReference with properties:

            XWorldLimits: [207000 209000]
            YWorldLimits: [911000 913000]
              RasterSize: [1000 1000]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 2
      CellExtentInWorldY: 2
    RasterExtentInWorldX: 2000
    RasterExtentInWorldY: 2000
        XIntrinsicLimits: [0.5 1000.5]
        YIntrinsicLimits: [0.5 1000.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []

Step 4: Visualize the Results

Display each image in map coordinates, and mark a reference point with a red + in both figures.

xlimits = [208000 208250];
ylimits = [911800 911950];

x = 208202.21;
y = 911862.70;

figure
mapshow(I_orig,R_orig)
hold on
plot(x,y,'r+')
xlim(xlimits)
ylim(ylimits)
ax = gca;
ax.TickDir = 'out';

figure
mapshow(I_half,R_half)
hold on
plot(x,y,'r+')
xlim(xlimits)
ylim(ylimits)
ax = gca;
ax.TickDir = 'out';

Graphically, they coincide, even though the same map location corresponds to two different locations in intrinsic coordinates.

[xIntrinsic1, yIntrinsic1] = worldToIntrinsic(R_orig, x, y)
xIntrinsic1 = 
       1202.7


yIntrinsic1 = 
       1137.8


[xIntrinsic2, yIntrinsic2] = worldToIntrinsic(R_half, x, y)
xIntrinsic2 = 
        601.6


yIntrinsic2 = 
       569.15


format(currentFormat);

Data Set Information

The concord_ortho_w.tif and concord_ortho_w.tfw files are derived using orthophoto tiles from the Bureau of Geographic Information (MassGIS), Commonwealth of Massachusetts, Executive Office of Technology and Security Services. For more information about the data sets, use the command type concord_ortho.txt. For an updated link to the data provided by MassGIS, see https://www.mass.gov/service-details/massgis-data-layers.

See Also

| | | | (Image Processing Toolbox)

Mapping Toolbox Documentation

Support