g_proj_GCP.m 2.57 KB
Newer Older
1
function [lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,frameRef)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%G_PROJ_GCP
%
% This function projects the GCPs that have an elevation greater than 0 
% onto the water level of zero elevation. This is a trick in order to be 
% able  to deal with GCPs that have elevation. The principle can be 
% illustrated in 1D like this. Consider the position x of a GCP that has an 
% elevation h. Then, the projection (xp) of this GCP onto the plane that has 
% zero elevation and along a line that joins the camera position of elevation
% H and the GCP is given by:
%
%      xp = x * (H - h)/H
%
% INPUTS:
%
16
17
%     LON0:     The longitiude or x-position of the camera
%     LAT0:     The latitude or y-position of the camera
18
%     H:        The elevation of the camera
19
20
%     lon_gcp0: The longitudes or x-position of the GCPs
%     lat_gcp0: The latitudes or y-position of the GCPs
21
22
23
%     h_gcp:    The elevations of the GCPs
%
% OUTPUTS:
24
25
%     lon_gcp:  The projected longitudes or x-position of the GCPs
%     lat_gcp:  The projected latitudes or y-position of the GCPs
26
27
28
29
%

n_gcp = length(h_gcp);

30
for i = 1:n_gcp    
31
    
32
    if strcmp(frameRef,'Geodetic') == true
33
34
35
36
37
38
39
40
41
42
43
        
        % The calculation for the distance between the GCPs and the camera
        % is done on an elliptical earth and this takes a long calculation
        % with the command m_idist from the m_map package. Therefore, the
        % calculation is done only on the GCPs with no elevation, hence
        % the 'if' condition here. 
        if h_gcp(i) > 0
            [range,a12,a21] = m_idist(LON0,LAT0,lon_gcp0(i),lat_gcp0(i));
            proj_factor     = H/(H-h_gcp(i));
            new_range       = range*proj_factor;
            [lon_gcp(i),lat_gcp(i),a21] = m_fdist(LON0,LAT0,a12,new_range);
44
45
46
            if lon_gcp(i) > 180
                lon_gcp(i) = lon_gcp(i) - 360;
            end
47
48
49
50
51
52
        else
            % If there's no elevation, then there's nothing to do.
            lon_gcp(i) = lon_gcp0(i);
            lat_gcp(i) = lat_gcp0(i);
        end
        
53
    elseif strcmp(frameRef,'Cartesian') == true
54
        
55
        % For cartesian situations where positions are provided in meters rather
56
57
58
59
60
61
62
63
64
65
66
67
68
        % than in latitude longitude the calculation is simpler.
        dx          = lon_gcp0(i) - LON0;
        dy          = lat_gcp0(i) - LAT0;
        range       = sqrt(dx^2 + dy^2);
        proj_factor = H/(H-h_gcp(i));
        new_range   = range*proj_factor;
        beta        = atan2(dy,dx);
        
        lon_gcp(i) = new_range*cos(beta) + LON0;
        lat_gcp(i) = new_range*sin(beta) + LAT0;
        
    end
end