How to drape Landsat images over Bedmap2 topography
This example combines a few built-in Matlab functions with some functions you'll find on the Mathworks File Exchange site.
To follow along with all the steps of this example, you'll need
Suppose you're interested in a 100-kilometer-wide area surrounding (77.1°S,162.7°E). Using the USGS lat/lon to path/row converter, we find that this area is covered by path 59, row 115. Let's initialize a polar stereographic map of our region and use the landsat function to get a quick look at the image we'll be working with. To plot the most recent Landsat image, you can simply use landat(59,115), but as I'm writing this example, the most recent relatively cloud-free image is from January 25, 2015. To ensure that you can get the same results I get, we'll specify landsat(59,115, 'Jan 25, 2015.'):
% Define region of interest: latcenter = -77.1; loncenter = 162.7; mapwidth = 100; % kilometers % Initialize a polarstereographic map: mapzoom(latcenter,loncenter,mapwidth) % Plot a low-res Landsat image of the area: landsat(59,115,'Jan 25, 2015')
That's a real nice photograph. Let's drape a high-resolution version of it over Bedmap2 topography. You'll need to log in to the Earth Explorer website and download the Level 1 GeoTIFF file.
After downloadting the GeoTiff, load the red, green, and blue bands individually with Aslak Grinsted's geoimread function. The GeoTIFF file contains georeferencing information in its metadata, so we can enter the center coordinates we've already defined as latcenter and loncenter. geoimread also accepts a "buffer" width, which corresponds to how many extra meters of data surrounding latcenter,loncenter we want to load. In this case, the buffer is equivalent to the half-width of our map, in meters. That is, buffer = mapwidth*1000/2.
% Load red band and polar stereographic coordinates of each pixel: [red,x,y] = geoimread('LC80591152015025LGN00_B4.TIF',... latcenter,loncenter,mapwidth*1000/2); % Load green band: green = geoimread('LC80591152015025LGN00_B3.TIF',... latcenter,loncenter,mapwidth*1000/2); % Load blue band: blue = geoimread('LC80591152015025LGN00_B2.TIF',... latcenter,loncenter,mapwidth*1000/2);
Concatenate red, green, and blue components into a single image:
rgb = cat(3,red,green,blue);
We can use bedmap2_interp to get the elevation of each pixel, but to do so, we must convert our x/y vectors into a gridded set of x/y coordinates, then convert gridded x/y into gridded lat/lon coordinates:
% Get gridded lat/lon corresponding to each image pixel: [xgrid,ygrid] = meshgrid(x,y); % Convert xgrid,ygrid to latgrid,longrid: [latgrid,longrid] = ps2ll(xgrid,ygrid); % Get surface elevation of each pixel: sfz = bedmap2_interp(latgrid,longrid,'surface','cubic');
The bedmap2_interp function returns NaN values over the ocean, so we need to set those NaNs to zero elevation:
To overlay an RGB image over a surface, you can use surf with the texturemap option, but that's somewhat unintuitive and can be a bit buggy. We'll use the Image Processing Toolbox's warp function instead. The only thing you have to remember with any Image Processing functions is, they typically plot in ij coordinates. That's why we have to switch back to xy coordinates with axis xy:
figure('pos',[100 100 1200 450]) warp(x,y,sfz,rgb) axis xy tight off view(-20,80);
The image above is okay, but shouldn't snow be white? (On a related note, Pro tip: if you ever do field work in snowy areas, disable the auto white balance on your camera, lest you you end up with gray snow). To enhance our Landsat image, use the Image Processing Toolbox functions imadjust and stretchlim:
rgbwow = imadjust(rgb,stretchlim(rgb),); figure('pos',[100 100 1200 450]) warp(x,y,sfz,rgbwow) axis xy tight off view(-20,80);
This example was written by Chad A. Greene of the University of Texas Institute for Geophysics (UTIG). February 2015.