readBasemapImage: How to transform the MapCellReference coordinates into non-projected coordinates

-- edit removed calculation to UTM and only display WGS-84 values

Dear Matlab Community,

I am currently working on loading tiles with ReadBaseMapImage and displaying them in a figure. However, I am facing an issue where the tiles are not aligned properly, and I'm unsure of the reason behind it.

I am restricted to using only a figure() for displaying the maps. I understand that there are other methods available in Matlab for displaying maps that are not part of the figure itself. However, my current approach follows these steps:

1. I have or create the corners of the map tiles in wgs84 coordinates.
2. I read the image and MapCellReference using [A,R,attrib] = readBasemapImage(basemap,latlim,lonlim)
3. I convert the #WorldLimits into WGS84 using [lat,lon] = projinv(proj,x,y).
4. I display both images where the world limits are located.

However, when I overlay these pictures, I notice that they are not properly aligned. It seems that there might be an error in my calculations.

My question is whether someone knows the source of this offset and if there is a way to fix it so that the tiles align with each other.

I would greatly appreciate any support or insights you can provide. Thank you in advance!

PS: The code I have used is as follows:

%% Getting Coordinates from a specific location in the world (e.g., Rome)
% Defining the points for the first tile
wgs_top_left_tile_1 = [41.891352, 12.490629];
wgs_bot_right_tile_1 = [41.888996, 12.494105];

% Defining the points for the second tile
wgs_top_left_tile_2 = [41.892278, 12.491798];
wgs_bot_right_tile_2 = [41.888096, 12.49414];

% Setting up the image type
map_type = 'satellite';

%% Calculation of the coordinates in WGS for tile 1
% Convert the input to limits
wgs_lim_lat_tile_1 = sort([wgs_top_left_tile_1(1), wgs_bot_right_tile_1(1)]);
wgs_lim_lng_tile_1 = sort([wgs_top_left_tile_1(2), wgs_bot_right_tile_1(2)]);

%% Calculation of the coordinates in WGS for tile 2
% Convert the input to limits
wgs_lim_lat_tile_2 = sort([wgs_top_left_tile_2(1), wgs_bot_right_tile_2(1)]);
wgs_lim_lng_tile_2 = sort([wgs_top_left_tile_2(2), wgs_bot_right_tile_2(2)]);

%% Get the map images
% Read the image for the first map
[my_image, img_info, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_1, wgs_lim_lng_tile_1);
% Read the Image for the second map
[my_image_2, img_info_2, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_2, wgs_lim_lng_tile_2);

%% Display the maps
% Create a figure
figure();

% Create axes for the figure
h = axes();

% Ensure Y axes are not inverted
h.YDir = 'normal';

% Enable hold on to draw multiple maps
hold on;

% Display the first map
display_image(my_image, img_info, h);

% Display the second map
display_image(my_image_2, img_info_2, h);

function display_image(img_rgb_matrix, img_map_info, my_axes)
    % Display an image with transparency on the specified axes, using WGS coordinates.

    % Create the image object
    x = image(img_rgb_matrix);

    % Set its transparency
    x.AlphaData = 0.5;

    % Convert pseudo-mercator coordinates to WGS coordinates
    [world_limits_x, world_limits_y] = get_wgs_from_psm(img_map_info.ProjectedCRS, img_map_info.XWorldLimits, img_map_info.YWorldLimits);

    % Set the x and y data for the picture
    x.XData = world_limits_y;
    x.YData = world_limits_x;

    % Set the x and y limits of the axes to the WGS world limits
    xlim(my_axes, world_limits_y);
    ylim(my_axes, world_limits_x);

    % Preserve the WGS scale by making the aspect ratio of the plot equal
    axis(my_axes, 'equal');
end

function [wgs_lat, wgs_lng] = get_wgs_from_psm(projection_crs, psm_x, psm_y)
    % Convert pseudo-mercator coordinates to WGS coordinates using the specified projection.

    % Convert projection coordinates to latitude and longitude
    [wgs_lat, wgs_lng] = projinv(projection_crs, psm_x, psm_y);
end

[Varun]

Hello! I found that manually adjusting the initial wgs coordinates can make the images overlap better. However, I am unsure about why this is the case. I used the coordinates below to get a clearer image.

% Defining the points for the first tile
wgs_top_left_tile_1 = [41.891352, 12.490629];
wgs_bot_right_tile_1 = [41.888996, 12.494105];

% Defining the points for the second tile
wgs_top_left_tile_2 = [41.892253, 12.491798];
wgs_bot_right_tile_2 = [41.888091, 12.49414];

Can you tell me more about how you arrived at your initial points for the tile? That might help me understand the issue better!

Answer 1

[Max]

Hi Varun,

thanks for your reply. The reason why your coordinates are woking is that the north/south middle of both pictures is the same vector.

The issue I had was that both pictures are provided flipped by readBasemapImage.

It is working by changing the code to:

%% Getting Coordinates from a specific location in the world (e.g., Rome)
% Defining the points for the first tile
wgs_top_left_tile_1 = [41.891352, 12.490629];
wgs_bot_right_tile_1 = [41.888996, 12.494105];
% Defining the points for the second tile
wgs_top_left_tile_2 = [41.892278, 12.491798];
wgs_bot_right_tile_2 = [41.888096, 12.49414];

% Setting up the image type
map_type = 'satellite';

%% Calculation of the coordinates in WGS for tile 1
% Convert the input to limits
wgs_lim_lat_tile_1 = sort([wgs_top_left_tile_1(1), wgs_bot_right_tile_1(1)]);
wgs_lim_lng_tile_1 = sort([wgs_top_left_tile_1(2), wgs_bot_right_tile_1(2)]);

%% Calculation of the coordinates in WGS for tile 2
% Convert the input to limits
wgs_lim_lat_tile_2 = sort([wgs_top_left_tile_2(1), wgs_bot_right_tile_2(1)]);
wgs_lim_lng_tile_2 = sort([wgs_top_left_tile_2(2), wgs_bot_right_tile_2(2)]);

%% Get the map images
% Read the image for the first map
[my_image, img_info, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_1, wgs_lim_lng_tile_1);
% Read the Image for the second map
[my_image_2, img_info_2, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_2, wgs_lim_lng_tile_2);

%% Display the maps
% Create a figure
figure();

% Create axes for the figure
h = axes();

% Ensure Y axes are not inverted
h.YDir = 'normal';

% Enable hold on to draw multiple maps
hold on;

% Display the first map
display_image(my_image, img_info, h);

% Display the second map
display_image(my_image_2, img_info_2, h);

function display_image(img_rgb_matrix, img_map_info, my_axes)
    % Display an image with transparency on the specified axes, using WGS coordinates.

    % Create the image object
    x = image(flipud(img_rgb_matrix));

    % Set its transparency
    x.AlphaData = 0.5;

    % Convert pseudo-mercator coordinates to WGS coordinates
    [world_limits_x, world_limits_y] = get_wgs_from_psm(img_map_info.ProjectedCRS, img_map_info.XWorldLimits, img_map_info.YWorldLimits);

    % Set the x and y data for the picture
    x.XData = world_limits_y;
    x.YData = world_limits_x;

    % Set the x and y limits of the axes to the WGS world limits
    xlim(my_axes, world_limits_y);
    ylim(my_axes, world_limits_x);

    % Preserve the WGS scale by making the aspect ratio of the plot equal
    axis(my_axes, 'equal');
end

function [wgs_lat, wgs_lng] = get_wgs_from_psm(projection_crs, psm_x, psm_y)
    % Convert pseudo-mercator coordinates to WGS coordinates using the specified projection.

    % Convert projection coordinates to latitude and longitude
    [wgs_lat, wgs_lng] = projinv(projection_crs, psm_x, psm_y);
end

Answer 2

[Max]

In addition, to controbute value for pepole with similar issues I created this class for transforming pictures also from Pseudo-Mercator to WGS or UTM. My initial idea was to transform these pictures to utm and this transformation is not as easy as the transformation to WGS. I hope this one helps at least one or two of you :)

These are the results:

The class can be executed with this script:

%% Getting Coordinates from a specific location in the world (e.g., Rome)
% Defining the points for the first tile
wgs_top_left_tile_1 = [41.891352, 12.490629];
wgs_bot_right_tile_1 = [41.888996, 12.494105];

% Defining the points for the second tile
wgs_top_left_tile_2 = [41.892278, 12.491798];
wgs_bot_right_tile_2 = [41.888096, 12.49414];

% Setting up the image type
map_type = 'satellite';

%% Calculation of the coordinates in WGS for tile 1
% Convert the input to limits
wgs_lim_lat_tile_1 = sort([wgs_top_left_tile_1(1), wgs_bot_right_tile_1(1)]);
wgs_lim_lng_tile_1 = sort([wgs_top_left_tile_1(2), wgs_bot_right_tile_1(2)]);

%% Calculation of the coordinates in WGS for tile 2
% Convert the input to limits
wgs_lim_lat_tile_2 = sort([wgs_top_left_tile_2(1), wgs_bot_right_tile_2(1)]);
wgs_lim_lng_tile_2 = sort([wgs_top_left_tile_2(2), wgs_bot_right_tile_2(2)]);

%% Get the map images
% Read the image for the first map
[my_image, img_info, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_1, wgs_lim_lng_tile_1);
% Read the Image for the second map
[my_image_2, img_info_2, ~] = readBasemapImage(map_type, wgs_lim_lat_tile_2, wgs_lim_lng_tile_2);

%% Display the maps
% Create a figure
fig = figure();

% Create axes for the figure
h = axes();

% Ensure Y axes are not inverted
h.YDir = 'normal';

% Enable hold on to draw multiple maps
hold on;

% Display the first map
display_image(my_image, img_info, h);

% Display the second map
display_image(my_image_2, img_info_2, h);

xlabel(h,'UTM easting') ;
ylabel(h,'UTM northing');
title(h,'ROME');

function display_image(img_rgb_matrix, img_map_info, my_axes)
    % Display an image with transparency on the specified axes, using WGS coordinates.
    image_transformer = class_maps_image_transformation();
    [wgs_pic, wgs_info, alpha] = image_transformer.transform_image(img_rgb_matrix, img_map_info, 'utm', 32);
    % Create the image object
    x = image(wgs_pic);

    % Set its transparency
    x.AlphaData = alpha*.5;

    


    % Set the x and y data for the picture
    x.XData = wgs_info.XWorldLimits;
    x.YData = wgs_info.YWorldLimits;

    % Set the x and y limits of the axes to the WGS world limits
    xlim(my_axes, wgs_info.XWorldLimits);
    ylim(my_axes, wgs_info.YWorldLimits);

    % Preserve the WGS scale by making the aspect ratio of the plot equal
    axis(my_axes, 'equal');
    
end

And this is the class used:

%--------------------------------------------------------------------------
% Title: MAPS IMAGE TRANSFORMATION
%--------------------------------------------------------------------------
% Description: A class for transforming an image from Web Mercator 
% projection to Universal Transverse Mercator coordinate system
%
% Compatibility: 
% The full functionality is only available on MATLAB R2020b or newer. 
% The following toolboxes are required:
% - 'Image Processing Toolbox'
% - 'Mapping Toolbox'
%
% License:
%    CC BY-SA 4.0
%    https://creativecommons.org/licenses/by-sa/4.0/
%
% Author: Max
% Date: 18.08.2023
%--------------------------------------------------------------------------
classdef class_maps_image_transformation
    %MAPS_IMAGE_TRANSFORMATION  Class for converting Web Mercator to UTM
    %pictures
    %
    % Properties (Access=protected, Constant)
    %    UTM_ZONE_ZERO - Constant for projcrs method
    %    FILL_VALUE_COLOR - Fill color for conversion which will be turned
    %                      to alpha afterwards
    %
    % Methods:
    %    class_maps_image_transformation - Construct an instance of this class
    %    transform_image - Transform the image of the tile from pseudo Mercator to UTM.
    %    check_for_toolboxes - Checks if all toolboxes are available.

    %------------ BEGIN OF PROPERTIES ------------
    properties (Access=protected, Constant)
        % Constant for projcrs method
        UTM_ZONE_ZERO int64 = 32600; 
        %Fill color for conversion which will be turned to alpha afterwards
        FILL_VALUE_COLOR = [255, 0, 0];
    end
    %------------- END OF PROPERTIES -------------
    %------------- BEGIN OF METHODS --------------
    methods
        function self = class_maps_image_transformation()
            %MAPS_IMAGE_TRANSFORMATION Construct an instance of this class.
            %    Check if all toolboxes are available
            self.check_for_toolboxes();
        end

        function [map_image, map_info, alpha_matrix] = ...
                transform_image(self, map_image_psm, map_info_psm, target_system, utm_zone)
            %TRANSFORM_IMAGE  Transform the image of the tile from Pseudo-Mercator to 
            %UTM.
            %This function applies a series of transformations to convert the tile's image 
            %from Pseudo-Mercator coordinates to UTM (Universal Transverse Mercator) 
            %coordinates, allowing it to be displayed correctly and aligned with the other
            %map data.
            %
            % Syntax:
            %    [new_img, new_info, alpha] = my_class.transform_image(image, info, 'utm', utm_zone);
            %    [new_img, new_info, alpha] = my_class.transform_image(image, info, 'wgs');   
            %
            % Inputs:
            %    self           - The class itself
            %    map_image_psm  - m*n*3 RGB image in the pseudo mercator
            %                     system
            %    map_info_psm   - map.rasterref.MapCellsReference
            %                     corresponding to the picture
            %    target_system  - System to convert to ('wgs' or 'utm')
            %    utm_zone       - only if targte system is 'utm'. Target utm zone
            %
            % Outputs:
            %    map_image_utm  - m*n*3 RGB image in the wgs system
            %    map_info_utm   - map.rasterref.MapCellsReference
            %                     corresponding to the wgs picture
            %    alpha_matrix   - m*n*1 matrix containing the picture alpha
            %
            % Example: 
            %    my_class = maps_image_transformation();
            %    [img, infos] = [my_image, img_info, ~] = readBasemapImage(map_type, lat, lng);
            %    [new_img, new_info, alpha] = my_class.transform_image(img, info, 'utm',34);
            %    [new_img, new_info, alpha] = my_class.transform_image(img, info, 'wgs');
            %    x = image(new_img);
            %    x.AlphaData = alpha*.5;
            %    x.XData = new_info.XWorldLimits;
            %    x.YData = new_info.YWorldLimits;
            %
            % Other m-files required: none
            % Subfunctions: none
            % MAT-files required: none
            %
            % See also: imref2d, projcrs, maprefcells, fitgeotrans,
            % fitgeotform2d, readBasemapImage
            
            %------------- BEGIN CODE --------------
            % Check if map_image_psm is a valid RGB image
            if isnumeric(map_image_psm) && ...
                    ndims(map_image_psm) == 3 &&...
                    size(map_image_psm, 3) == 3
                % Do nothing
            else
                error('map_image_psm is not a valid m-by-n-by-3 RGB matrix.');
            end
            % Check if map_info_psm is of type map.rasterref.MapCellsReference
            if isa(map_info_psm, 'map.rasterref.MapCellsReference')
                % Do nothing
            else
                error('map_info_psm is not of type map.rasterref.MapCellsReference.');
            end

            if ischar(target_system) && (strcmp(target_system, 'utm') || strcmp(target_system, 'wgs'))
                % Do Nothing
                
                if strcmp(target_system, 'utm')
                    % Check if UTM zone is provided
                    if nargin < 5
                        error('UTM zone is required for target_system = ''utm''.');
                    end
                    
                    % Additional checks for UTM zone, you can modify this part as needed
                    if ~isnumeric(utm_zone) || ~isscalar(utm_zone) || utm_zone < 1 || utm_zone > 60
                        error('Invalid UTM zone. UTM zone should be a numeric scalar between 1 and 60.');
                    end
                    
                elseif strcmp(target_system, 'wgs')
                    % Do Nothing
                end
                
                % Here you can continue with the rest of your transformation logic
                
            else
                error('target_system is not valid. It should be either ''utm'' or ''wgs''.');
            end
            
            % Read out the pseudo Mercator transformation points
            psm_x_lims = map_info_psm.XWorldLimits;
            psm_y_lims = map_info_psm.YWorldLimits;
            
            % Generates a box around the picture in pseudo mercator
            psm_points_x = [ ...
                min(psm_x_lims), ...
                min(psm_x_lims), ...
                max(psm_x_lims), ...
                max(psm_x_lims) ...
                ];
            psm_points_y = [ ...
                min(psm_y_lims), ...
                max(psm_y_lims), ...
                max(psm_y_lims), ...
                min(psm_y_lims), ...
                ];
            
            % Transform this box from PSM to WGS84
            [wgs_lat_lims, wgs_lng_lims] = projinv( ...
                            map_info_psm.ProjectedCRS, ...
                            psm_points_x, ...
                            psm_points_y...
                            );
            if(strcmp(target_system, 'utm'))
                % Transform this box from WGS zu UTM
                if verLessThan('MATLAB', '9.9') % before MATLAB R2020b (Version 9.9)
                    error('The map base plugin requires MATLAB R2020b or newer');
                end
                target_crs = projcrs(self.UTM_ZONE_ZERO + utm_zone);
                % take the utm values for final values
                [utm_easting, utm_northing] = projfwd(...
                                target_crs, ...
                                wgs_lat_lims, ...
                                wgs_lng_lims...
                                );
                % Control points 
                control_target = transpose([utm_easting; utm_northing]);
            
                % Control points PSM
                control_psm = transpose([psm_points_x; psm_points_y]);
                
                % Generate the image reference for the not transfored image
                RA = imref2d(...
                    map_info_psm.RasterSize, ...
                    map_info_psm.XWorldLimits, ...
                    map_info_psm.YWorldLimits);
                
                % Flip the image
                flipped_image = flipud(map_image_psm);
            
                % Generate a transforming class
                if verLessThan('MATLAB', '9.13') % before MATLAB R2022b (Version 9.13)
                    % use fitgeotrans (slower)
                    t_raster2wgs = fitgeotrans(control_psm, control_target, 'similarity');
                else
                    % use fitgeotform2d (faster)
                    t_raster2wgs = fitgeotform2d(control_psm, control_target, 'similarity');
                end
                
                % Transform the image
                [utm_image, RB]= imwarp(flipped_image, RA, t_raster2wgs, 'FillValues', ...
                    self.FILL_VALUE_COLOR);
            
                % Create the image object
                map_image = utm_image;
                
                % store the x and y data of the utm image
                tile_x_data = [min(RB.XWorldLimits), max(RB.XWorldLimits)];
                tile_y_data = [min(RB.YWorldLimits), max(RB.YWorldLimits)];
    
                % store image information
                map_info = maprefcells(...
                    tile_x_data, ...
                    tile_y_data, ...
                    size(utm_image));
                % Create new alpha matrix
                alpha_matrix = ones([size(map_image,1), size(map_image,2)]) ;
    
                % set alpha matrix to transparent were the Fill color is met
                alpha_matrix((map_image(:,:,1) == self.FILL_VALUE_COLOR(1)) & ...
                          (map_image(:,:,2) == self.FILL_VALUE_COLOR(2)) & ...
                          (map_image(:,:,3) == self.FILL_VALUE_COLOR(3))) = 0;
            else
                tile_x_data = [min(wgs_lng_lims), max(wgs_lng_lims)];
                tile_y_data = [min(wgs_lat_lims), max(wgs_lat_lims)];
                map_info = maprefcells(...
                    tile_x_data, ...
                    tile_y_data, ...
                    size(map_image_psm));
                % Create the image object
                map_image = flipud(map_image_psm);
                % Create new alpha matrix
                alpha_matrix = ones([size(map_image,1), size(map_image,2)]) ;
            end
        
            
            %------------- END OF CODE --------------
        end
        
        function check_for_toolboxes(~)
            %CHECK_FOR_TOOLBOXES  Check if all needed toolboxes are
            %available
			%Check if the mapping and image toolbox are available on
            %MATLAB. Throws and error if toolboxes are not available.
			%
			% Inputs:
			%    None.
			%
			% Outputs:
			%    None.
			%
			% Example:
			%    self.check_for_toolboxes()
            %
            % See also: matlab.addons.installedAddons, strcmp
            
            %------------- BEGIN CODE --------------
            %% Checks for the image and mapping toolbox
            % Name of the toolbox you want to check if it's installed
            toolboxName_image = 'Image Processing Toolbox';
            toolboxName_mapping = 'Mapping Toolbox';
            
            % Get information about all installed add-ons
            installedAddons = matlab.addons.installedAddons();
            
            % Initialize a variable to check if the toolbox is installed
            isToolboxInstalled_image    = false;
            isToolboxInstalled_mapping  = false;
            
            % Loop through the installed add-ons to find the toolbox
            for i = 1:size(installedAddons, 1)
                if strcmp(installedAddons.Name(i), toolboxName_image)
                    isToolboxInstalled_image = true;
                end
                if strcmp(installedAddons.Name(i), toolboxName_mapping)
                    isToolboxInstalled_mapping = true;
                end
                if (isToolboxInstalled_mapping == true) && ...
                        (isToolboxInstalled_image == true)
                    break;
                end
            end
            if (isToolboxInstalled_image == false)
                error('The image toolbox is required for image transformation');
            end
            if (isToolboxInstalled_mapping == false)
                error('The mapping toolbox is required for image transformation');
            end
            %------------- END OF CODE --------------
        end
    end
    %------------- END OF METHODS --------------
end

[Juan Sebastian]

Thanks a lot!!