Geographic Information System (GIS) Using MATLAB
MATLAB provides robust support for Geographic Information System (GIS) operations. With specialized toolboxes such as the Mapping Toolbox, users can process geospatial data, visualize maps, and solve geospatial problems effectively.
Introduction to GIS in MATLAB
GIS involves the collection, processing, and analysis of geospatial data. MATLAB supports various file formats, including shapefiles (.shp), GeoTIFF (.tif), and KML (.kml).
Working with Shapefiles
Shapefiles are one of the most commonly used GIS data formats. MATLAB provides the shaperead function to load shapefiles.
% Reading a shapefile
S = shaperead('countries.shp');
disp(S); % Display shapefile structure
To visualize the shapefile, use the mapshow function:
% Visualizing shapefile data
figure;
mapshow(S);
title('World Map');
Output: Displays a map of countries with their boundaries.
Geographic Coordinates and Projections
GIS data often uses geographic coordinates (latitude and longitude). MATLAB provides functions for converting between different coordinate systems.
Example: Converting Geographic Coordinates to Cartesian Coordinates
% Convert geographic to Cartesian coordinates
lat = 37.7749; % Latitude of San Francisco
lon = -122.4194; % Longitude of San Francisco
[East, North, Up] = geodetic2enu(lat, lon, 0, 0, 0, 0, wgs84Ellipsoid);
disp([East, North, Up]);
Output: Displays the Cartesian coordinates of the input location.
Plotting Geographic Data
MATLAB provides various tools to visualize geospatial data on maps.
Example: Creating a Geographic Plot
% Plotting geographic data
latitudes = [37.7749, 34.0522, 40.7128]; % San Francisco, LA, NY
longitudes = [-122.4194, -118.2437, -74.0060];
geoplot(latitudes, longitudes, 'r*-');
geobasemap('streets');
title('Major US Cities');
Output: A geographic plot of the given locations with a street map background.
Raster Data Processing
Raster data represents information in grid format (e.g., elevation or temperature data). MATLAB supports raster processing using the readgeoraster function.
Example: Reading and Displaying Raster Data
% Reading raster data
[Z, R] = readgeoraster('elevation.tif');
disp(Z); % Elevation data grid
Output: Displays elevation values as a grid.
To display the raster data:
% Visualizing raster data
figure;
mapshow(Z, R, 'DisplayType', 'surface');
colormap('jet');
colorbar;
title('Elevation Map');
Output: A color-coded map representing elevation data.
Advanced GIS Example Problems and Solutions in MATLAB
In this section, we delve into advanced GIS operations, including spatial analysis, map projections, geospatial statistics, and more.
Example 1: Spatial Buffering
Spatial buffering is a technique to create zones around geographic features. MATLAB provides functions to perform such operations.
% Example: Create a buffer around points
S = shaperead('cities.shp'); % Read shapefile
lat = [S.Lat]; % Extract latitudes
lon = [S.Lon]; % Extract longitudes
bufferRadius = 0.5; % Buffer radius in degrees
[latBuff, lonBuff] = bufferm(lat, lon, bufferRadius); % Create buffer
figure;
geoshow(lat, lon, 'DisplayType', 'point');
hold on;
geoshow(latBuff, lonBuff, 'DisplayType', 'polygon', 'FaceAlpha', 0.3);
title('Spatial Buffer Around Cities');
Output: A map showing points for cities and translucent buffer zones around them.
Example 2: Raster Algebra
Raster algebra involves mathematical operations on raster datasets. This is useful in applications like calculating slope or merging multiple layers.
% Example: Calculate slope from elevation data
[Z, R] = readgeoraster('elevation.tif');
[dZdx, dZdy] = gradient(Z); % Compute gradients
slope = atand(sqrt(dZdx.^2 + dZdy.^2)); % Calculate slope in degrees
figure;
imagesc(slope); % Visualize slope
colormap('hot');
colorbar;
title('Slope Map');
Output: A heat map representing the slope across the elevation data.
Example 3: Map Projection and Reprojection
Map projections are used to convert 3D geographic coordinates into 2D map representations. MATLAB supports multiple projection systems.
% Example: Reproject coordinates to UTM
lat = 37.7749; % Latitude
lon = -122.4194; % Longitude
utmZone = '10T'; % UTM zone
[x, y] = projfwd(projcrs(utmZone), lat, lon); % Convert to UTM
disp([x, y]);
Output: Displays UTM (x, y) coordinates for the given location.
Example 4: Interpolation of Geographic Data
Interpolation is often used to estimate values at unmeasured locations. MATLAB provides powerful interpolation tools for geospatial data.
% Example: Interpolate temperature data
lat = [34, 35, 36, 37]; % Latitudes of stations
lon = [-120, -121, -122, -123]; % Longitudes of stations
temp = [15, 18, 20, 17]; % Temperature at stations
[gridLat, gridLon] = meshgrid(34:0.1:37, -123:0.1:-120); % Grid for interpolation
gridTemp = griddata(lat, lon, temp, gridLat, gridLon, 'cubic'); % Interpolation
figure;
geoshow(gridLat, gridLon, gridTemp, 'DisplayType', 'surface');
colorbar;
title('Interpolated Temperature Map');
Output: A surface map displaying interpolated temperature values.
Example 5: Spatial Statistics
Performing spatial statistics helps analyze patterns and distributions of geospatial data.
% Example: Compute centroid of a polygon
S = shaperead('state_boundaries.shp'); % Read shapefile
[latC, lonC] = polygeom(S(1).X, S(1).Y); % Calculate centroid
disp([latC, lonC]);
Output: Displays the centroid coordinates of the first polygon in the shapefile.
Example 6: Route Optimization
Route optimization is a critical GIS application in transportation and logistics. MATLAB can compute shortest paths using graph theory.
% Example: Shortest path between two points
lat = [37.7749, 36.7783, 34.0522]; % San Francisco, Fresno, Los Angeles
lon = [-122.4194, -119.4179, -118.2437];
distances = [0, 3.5, 6.2; 3.5, 0, 2.8; 6.2, 2.8, 0]; % Distance matrix
G = graph(distances);
shortestPath = shortestpath(G, 1, 3); % Find shortest path from SF to LA
disp(shortestPath);
Output: Displays the shortest path through nodes [1 2 3].
Useful MATLAB Functions for GIS
Practice Questions
Test Yourself
1. Read a shapefile and plot its data using mapshow.
2. Convert the geographic coordinates of Los Angeles to Cartesian coordinates using geodetic2enu.
3. Read a GeoTIFF file containing temperature data and visualize it as a surface plot.
4. Create a geographic plot of three cities and display them on a satellite basemap.
5. Create a buffer zone around the borders of a country and visualize it on a map.
6. Compute the slope and aspect of a given elevation dataset.
7. Reproject coordinates from WGS84 to Mercator projection.
8. Interpolate precipitation data across a grid and visualize it as a contour plot.
9. Determine the shortest route between three major cities using a distance matrix.