Change Projection of Basemap Image

Change Projection Using Projection Stored in axesm-Based Map

This example shows how to project a basemap image by using the projection method stored in an axesm-based map. For more information about projections supported by axesm-based maps, see Summary and Guide to Projections.

Load Data

Read the latitude and longitude coordinates of European capitals from a text file.

[lat,lon] = readvars("european_capitals.txt");

Read Basemap Image

Specify latitude and longitude limits for the basemap image. Find the limits for the region containing the capitals by using the geoquadline function. Expand the limits by 10 degrees in latitude and longitude by using the bufgeoquad function.

[latlim,lonlim] = geoquadline(lat,lon);
[latlim,lonlim] = bufgeoquad(latlim,lonlim,10,10);

Read an image for the region from the "bluegreen" basemap, using a zoom level of 3, as an array and a map cells reference object in Web Mercator coordinates.

[A,RA] = readBasemapImage("bluegreen",latlim,lonlim,3);

Unproject Basemap Image

Extract the Web Mercator coordinates of the basemap image from the reference object. Then, unproject the coordinates to geographic coordinates.

[xGrid,yGrid] = worldGrid(RA);
[latGrid,lonGrid] = projinv(RA.ProjectedCRS,xGrid,yGrid);

The coordinates of the European capitals are already in geographic coordinates. If you are adapting this example to your own code and your data is in projected coordinates, you must unproject your data to geographic coordinates.

Create Map

Create an axesm-based map with a projection that is suitable for the region by using the worldmap function. Alternatively, you can create an axesm-based map with a specified projection ID by using the axesm function.

figure
h = worldmap(latlim,lonlim);

View the projection method stored in the map.

getm(h,"mapprojection")

ans = 
'eqdconic'

The result 'eqdconic' means that the map uses an Equidistant Conic projection.

Display the basemap and the European capitals on the map. The geoshow function projects and displays the latitude-longitude coordinates of the basemap image and capitals using the projection method stored in the map.

geoshow(latGrid,lonGrid,A)
geoshow(lat,lon,DisplayType="point",Marker="pentagram")
title("European Capitals")
subtitle("Basemap Attribution: Natural Earth")

Text in the projected basemap image can appear distorted.

Change Projection Using Projected CRS Object

This example shows how to read a basemap image using a projected CRS object and then display the projected image on an axes.

Load Data

Read a shapefile containing the coordinates of world cities into the workspace as a geospatial table. Extract the latitude and longitude coordinates.

GT = readgeotable("worldcities.shp");
lat = GT.Shape.Latitude;
lon = GT.Shape.Longitude;

Specify Projected CRS and Find Area of Use

Specify the projected CRS you want to use. For this example, create a projcrs object for WGS 84 / Equal Earth Asia-Pacific, which has the EPSG code 8859.

equalEarth = projcrs(8859);

Different projected CRSs are valid for different areas of use. For many projected CRSs, you can find the area of use by displaying the well-known text (WKT) string of the CRS and finding the BBOX attribute. This WKT specifies BBOX in the second-to-last line.

wktstring(equalEarth,Format="formatted")

ans = 
    "PROJCRS["WGS 84 / Equal Earth Asia-Pacific",
         BASEGEOGCRS["WGS 84",
             DATUM["World Geodetic System 1984",
                 ELLIPSOID["WGS 84",6378137,298.257223563,
                     LENGTHUNIT["metre",1]]],
             PRIMEM["Greenwich",0,
                 ANGLEUNIT["degree",0.0174532925199433]],
             ID["EPSG",4326]],
         CONVERSION["Equal Earth Asia-Pacific",
             METHOD["Equal Earth",
                 ID["EPSG",1078]],
             PARAMETER["Longitude of natural origin",150,
                 ANGLEUNIT["degree",0.0174532925199433],
                 ID["EPSG",8802]],
             PARAMETER["False easting",0,
                 LENGTHUNIT["metre",1],
                 ID["EPSG",8806]],
             PARAMETER["False northing",0,
                 LENGTHUNIT["metre",1],
                 ID["EPSG",8807]]],
         CS[Cartesian,2],
             AXIS["(E)",east,
                 ORDER[1],
                 LENGTHUNIT["metre",1]],
             AXIS["(N)",north,
                 ORDER[2],
                 LENGTHUNIT["metre",1]],
         USAGE[
             SCOPE["Very small scale equal-area mapping - Asia-Pacific-centred."],
             AREA["World centred on Asia-Pacific."],
             BBOX[-90,-29.99,90,-30.01]],
         ID["EPSG",8859]]"

Specify the latitude and longitude limits for the basemap image by using the information in BBOX. For this projected CRS, unwrap the maximum longitude by adding 360 degrees to the minimum longitude.

latlim = [-90 90];
lonmin = -29.99;
lonlim = [lonmin lonmin+360];

Read Basemap Image

Read an image from the "satellite" basemap as an array, a map cells reference object in Web Mercator coordinates, and an attribution string.

[A,RA,attrib] = readBasemapImage("satellite",latlim,lonlim);

Reproject Data and Basemap Image

Reproject the coordinates of the basemap image to Equal Earth coordinates.

[xGrid,yGrid] = worldGrid(RA);
[latGrid,lonGrid] = projinv(RA.ProjectedCRS,xGrid,yGrid);
[xEqualEarth,yEqualEarth] = projfwd(equalEarth,latGrid,lonGrid);

Project the geographic coordinates of the world cities to Equal Earth coordinates. If you are adapting this example to your own code and your data is in projected coordinates, you must first unproject your data to geographic coordinates.

[xData,yData] = projfwd(equalEarth,lat,lon);

Create Map

Display the projected basemap image and world cities on a map with no axis labels.

figure
mapshow(xEqualEarth,yEqualEarth,A)
hold on
mapshow(xData,yData,DisplayType="point",Marker=".")
title("World Cities")
subtitle("Basemap Attribution: " + attrib)
axis off

Text in the projected basemap image can appear distorted.