projinv - Inverse map projection using PROJ.4 map projection library

Syntax

Description

[lat, lon] = projinv(proj, x, y) returns the latitude and longitude values from the inverse projection transformation. proj is a structure defining the map projection. proj can be a map projection mstruct or a GeoTIFF info structure. x and y are x-y map coordinate arrays. For a complete list of GeoTIFF info and map projection structures that you can use with projinv, see the reference page for projlist.

Examples

Display Boston Orthophoto on a Mercator projection

1. Import the Boston roads from the shapefile and obtain the projection structure from the 'boston.tif' orthophoto:

   roads = shaperead('boston_roads.shp');
   proj = geotiffinfo('boston.tif');
   

2. Convert the road coordinates to the projection's length unit. As shown by the UOMLength field of the projection structure, the units of length in the projected coordinate system is US Survey Feet. Coordinates in the roads shapefile are in meters:

   proj.UOMLength
   
   ans =
   US survey foot
   
   x = [roads.X] * unitsratio('survey feet','meter');
   y = [roads.Y] * unitsratio('survey feet','meter');
   
   % Now convert the scaled coordinates of the roads 
   % to latitude and longitude.
   [roadsLat, roadsLon] = projinv(proj, x, y);
   

3. Read the boston_ovr.jpg image and worldfile:

   RGB = imread('boston_ovr.jpg');
   R = worldfileread(getworldfilename('boston_ovr.jpg'));

4. Read state boundary vectors for Massachusetts from the usastatehi shapefile using a selector to eliminate other states:

   S = shaperead('usastatehi', 'UseGeoCoords', true, ...
       'Selector',{@(name) strcmpi(name,'Massachusetts'), 'Name'});
   

5. Open a figure with a Mercator projection and display the state boundary, image, and roads:

   figure
   axesm('mercator')
   
   geoshow(S.Lat, S.Lon, 'Color','red')
   geoshow(RGB, R)
   geoshow(roadsLat, roadsLon, 'Color', 'green')
   

6. Set the map boundary to the image's northern, western, and southern limits, and the eastern limit of the state boandary within the image latitude bounding box:

   [lon, lat] = mapoutline(R, size(RGB(:,:,1)));
   ltvals = find((S.Lat>=min(lat(:))) & (S.Lat<=max(lat(:))));
   setm(gca,'maplonlimit',[min(lon(:)) max(S.Lon(ltvals))], ...
            'maplatlimit',[min(lat(:)) max(lat(:))])
   tightmap

   

   boston_ovr.jpg image copyright © GeoEye, all rights reserved.

See Also

geotiffinfo | mfwdtran | minvtran | projfwd | projlist

© 1984-2012- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS