Visualizing Raster Data in 3D using Matlab

Prepared By Joshua Roundy

Project for:

Geographical Information Systems in Water Recourses

5 Dec 2008

Report Prepared For: Dr. David G. Tarboton, Dr. David R. Maidment, and Dr. Ayse Irmak

Introduction

Viewing a time series of raster data allows for visualizing how a variable changes over time and space. Viewing the data in two dimensions allows for sufficient analysis of the x and y spatial change, however it does not allow for the change in elevation. Being able to view raster data with respect to three dimensions allows further analysis of how the variables change with elevation. The purpose of this project is to explore the basic functions of the Matlab Mapping Toolbox and develop Matlab scripts that can be used to view and analyze time series of raster data in 3 dimensions. These scripts are then used to perform some basic analysis of radar precipitation data over the Animas River Basin located near Durango, Colorado. Matlab was chosen as the development environment because of its strong 3D mapping capabilities and its effectiveness in reading in large datasets and the ability to create AVI movies. Another advantage of using Matlab that will not be explored in this project is its ability to directly link with complex models that are written in C or FORTRAN. The Matlab mapping toolbox was used extensively throughout this project, as such a quick summary of the toolbox will be provided.

Matlab Mapping Toolbox Overview

The Matlab Mapping toolbox is a compilation of functions used to work with mapping data.

It is an additional feature that does not come with the standard version of Matlab and musts be purchased separately. Detailed documentation about the toolbox and all the functions available is at the MathWorks website (MathWorks, 2008). Matlab has a specified format for storing both vector and raster data. The vector data is stored in a map structure that details all the information associated to the data. An example of this is shown in figure 1.

Figure 1- Matlab Map Structure

As can be seen from figure 1, the shown data is an example of polygon data named Area. Each polygon in the structure is stored by the index of the Area structure. For example the first polygon stored in the structure Area can be accessed by typing “Area (1)”. This holds true for all other vector data sets. The associated latitudes and longitudes for the polygon “Area (1)” are stored in the variables X and Y (Note the X and Y values are the x and y coordinates if the data is projected). These values can be accessed by typing “Area (1).X”. It is not necessary to use this map structure when creating maps as any vectors of x and y data can be mapped; however, it is an organized way of using map data and has the added advantage of being able to use other built in Mapping Toolbox functions. The simplest way to access vector data in Matlab is to use the shaperead function. In its simplest form the function takes only the name of the shapefile and returns a map structure of the data. One limitation to the shaperead function that was discovered in this project is its inability to read three dimensional vector data sets. To read such data sets they must be converted to normal polygon, polyline and point files. The function shapewrite does just the opposite of shaperead, as it takes a map structure and writes a shapefile. These simple functions provide and easy way of using vector data in Matlab.

It is also straight forward to use grid data in Matlab. There are several different ways of storing grid data in Matlab, however only the easiest and most common is used in this project. The grid data is stored in a normal Matlab matrix, but in order for it to be used in mapping it also must have an accompanying reference vector. The matrix together with the reference vector makes up the structure for storing raster data in Matlab. The reference vector contains three elements, the number of cells per degree, the north latitude, and the west longitude of the grid. The first row in the matrix is the southernmost row of the data.

Figure 2 - Raster Data in Matlab

Figure 2 illustrates the reference vector and matrix relationship, where the red represents the location of the reference vector and its corresponding grid in the matrix. When working with raster data in Matlab it is important to remember this relationship so the grid data is not viewed or analyzed upside down. There are built-in Matlab functions for reading in grid data in different formats, some of which will be covered in detail later on, however it is not mandatory that these functions are used. A valid raster data set can be any matrix that is read or created in Matlab as long as it has a reference vector.

Once the data is in Matlab there are several ways to plot it. When plotting data in Matlab it is easiest to keep all data in latitude and longitude coordinates. The projection that the data is viewed in is defined when creating the axes, thus the data is not actually projected. There is a way of actually projecting the data; however this was not used in this project. Once a figure and axis are created with the appropriate projection, then there are specific functions used to plot the data. The easiest and most versatile function for plotting data is geoshow. The function geoshow can be used with both a vector and raster data and also takes a projection argument. Another helpful mapping function is meshm, which is used for grid data and allows the ability to overlay other datasets.

The Matlab mapping toolbox has hundreds of useful functions for analyzing GIS data. To this point only the most simple and basics functions were discussed to enable one to start using Matlab. For more details the reader is encouraged to explore the MathWorks website where there is more detailed and useful examples (MathWorks, 2008).

Developing Matlab Scripts

The purposes of this project were not only to learn about the Matlab Mapping Toolbox, but to write scripts in Matlab to analyze radar precipitation data with respect to elevation over the Animas River Basin, near Durango, Colorado. The Animas River Basin is located in both Colorado and New Mexico as shown in figure 3.

Figure 3 - Animas River Basin

The Animas River Basin is approximately 3,555 square kilometers and the Animas River has an average annual runoff for the last four years of 890 CFS (USGS Real-Time Water Data for the Nation). As can be seen from the figure, the basin also provides an interesting study of precipitation and elevation because it has a mountainous area along with a valley area.

The two types of data needed for this project were Digital Elevation Model (DEM) and precipitation radar data. The DEM data was obtained from the USGS website; specifically the one degree quadrangle in DEM format was download (USGS Geographic Data Download). This DEM data has a grid size of approximately 90 meters. This data was chosen because of the built in Matlab function that reads the data. The Matlab function usgsdems gives the names of the DEM files needed for a given latitude-longitude bounding box, which makes finding the needed files really simple. Once the appropriate files are downloaded, they are just as easily read into Matalb with the usgsdem function. This function returns a matrix of the data and reference vector. While working with this data it was discovered that for plotting purposes it was really slow if the data was read in at full resolution. The usgsdem function also takes a value for what resolution to read the data in at. The Animas River Basin covers six different quadrangles, as such a function to combine the six grids into one was written. Overall it was very simple to work with the DEM data because of the built in Matlab functions.

The radar precipitation data was obtained from the National Weather Service and was supplied by my major professor. The NEXRAD data is in XMRG binary format in the HRAP projection with six hour precipitation data for 1996 and 1997. The original goal of the project was to read the data in from the binary files, however after many attempts to do this failed, the data was obtained in ASCII format instead. A Matlab function was written to read the ASCII files into Matlab. Once the grid data was read into Matlab the problem of the HRAP projection was considered.

The HRAP projection uses a spherical datum with a Polar Stereographic projection that is shifted so that all grids have positive values in the United States. In order to properly transform the precipitation data to latitude and longitude coordinates the precipitation data was resampled. This was done by creating a grid in latitude and longitude coordinates with a grid size of approximately 2 kilometers. The centers of each grid were determined and the coordinates were converted to the HRAP projection and the radar data was sampled accordingly. In order for this procedure to work a function to convert from latitude and longitude to the HRAP projection was written. This was accomplished by following the outline given in the paper by Reed and Maidment which was very helpful (1999).

Once this procedure for reading the data was setup a routine was written to read in the data in at every six-hour period for a given time frame and store all grids into one Matlab Structure. This structure provided a way to compile all the grid data over a time period into one variable making it easier to perform statistical analysis and visualize the data. The statistical analysis performed was to calculate the total days of a wet grid, the average precipitation for the wet grid and the total precipitation over a certain time period. The function that was written to do this looped through the radar data structure computing each of the desired statistics. This was very easy to do as it used simple logical queries, matrix addition and scalar multiplication applied to the created structure.

A function for creating a movie of the precipitation data was also written as part of this project. The first step for doing this was to map the precipitation data over the DEM data. Using the meshm function the DEM data was mapped and scaled appropriately. Then the precipitation grid was overlaid on the DEM data. An example of this is shown in figure 4.

Figure 4 - Example of Precipitation Data overlaid on DEM data

This procedure was then recursively executed for each time step and saves each image into an AVI movie. The Matlab function for making movies is not unique to the Mapping Toolbox; as such it can be used by a normal version of Matlab. When the Matlab function gets the image for the movie it does it essentially like a print screen, so if any other object opens on top of the Matlab figure while running, it will also be saved in the image for the movie.

Results

The precipitation over the Animas River Basin was evaluated for the month of January for the year 1997.

Figure 5 - Total Wet Grids for January 1997

Figure 5 shows the total number of times the grid cell was wet during this time period. As can be seen from this figure the higher elevations have more wet grids than the valleys. The east side of the mountain range also seems to be wetter than the west side. Figure 6 shows the average amount of precipitation for the wet grids.

Figure 6 - Average Precipitation for Wet Grids for January 1997

As can be seen from the figure there are certain areas that receive more precipitation on average. The interesting thing here is that in comparison to wet grids there is a much smaller area of high precipitation. This shows that even though the higher elevations had similar number of times of wet grid the average amount of precipitation is very different. The area of dark blue in figure represents the part of the basin that receives the most precipitation during a wet period. Figure 7 shows the total precipitation over the basin.

Figure 7 - Total Precipitation for January 1997

From this figure it can be seen that the area of highest average precipitation does receive the most precipitation in the basin. Overall there appears to be a trend of increased precipitation with elevation.

To better examine the relationship of precipitation with elevation, the elevation of each grid was sampled by nearest neighbor of the center of the grid from the DEM data. A plot of the elevation vs. total precipitation is plotted in figure 8.

Figure 8 - Elevation vs. Precipitation

This shows that there is a definite trend of precipitation to elevation. This is particularly important for estimating the precipitation for a rainfall runoff model over a mountainous area. If a simple interpolation from rain gages is used for a mountainous area, then it is likely that the precipitation would be over or under estimated depending on the location of the rain gages. Using radar data is not without its share of uncertainties. The precipitation data analyzed in this study most likely fell as snow since it occurred in the month of January, so there is a question of the accuracy of the radar to properly measure snow. Also, because of the mountainous area the radar is measured far above the land surface which questions the accuracy of the precipitation data that actually fell to the ground. The best approach would be to combine the radar data with that of the rain gages in the basin. This would provide a way to use the rain gages to test and adjust the accuracy of the radar data which more correctly distributes the precipitation over the basin.

Conclusion

This project successfully showed the basic functions used in the Matlab Mapping Toolbox and utilized these and other functions to develop scripts to overlay precipitation data on a DEM. Matlab proved to be an excellent tool for creating useful maps and movies for visualizing this data. These tools were then used to analyze the relationship of elevation to precipitation over the Animas River Basin. The analysis showed that there is an increase in precipitation with an increase of elevation. There is also a lot of uncertainty to the accuracy of the radar data. To help with this accuracy, it should be combined with rain gages.

References

MathWorks (Accessed September 25), Mapping Toolbox 2.7.1

http://www.mathworks.com/products/mapping/

Reed, S.M., and D.R. Maidment, "Coordinate Transformations for Using NEXRAD Data in GIS-based

Hydrologic Modeling," Journal of Hydrologic Engineering, 4, 2, 174-182, April 1999

USGS (United States Geological Survey) (Accessed 2008, November 11), USGS Geographic Data

Download. http://edc2.usgs.gov/geodata/index.php

USGS (United States Geological Survey) (Accessed 2008, November 24), USGS Real-Time Water Data for

the Nation. http://waterdata.usgs.gov/nwis/rt