| Mapping Toolbox | |
Importing Other Data
With the huge quantities of geographic information available, it is impossible to provide interface functions for every file format. To help you read data in other formats, the Mapping Toolbox provides a set of high-level file import functions. These allow you to read data from binary or text files by simply describing the contents of the file. The files may contain matrices or fixed length records of geographic data. MATLAB also provides some file import functions. Type help iofun for more information on the MATLAB functions.
Much geographic data is provided as very large files containing matrices of values. You can read the part of the matrix of interest to you using the readmtx function.
To illustrate the use of readmtx, read in data giving the intensity of stable light sources collected from the Defense Meteorological Satellite program <http://www.ngdc.noaa.gov/dmsp/dmsp.html>. Download a recent data file (e.g., usa.bil) and header file (e.g., usa.hdr) from the DMSP FTP site. We will read and display data for the Cape Cod region. Take the following minimum information from the header file: the number of rows (NROWS) and columns (NCOLS) in the matrix, the location of the upper left corner (ULXMAP, ULYMAP), the grid spacing (XDIM, YDIM), and the storage format of the numbers (NBITS).
Most matrix files can be read using the same procedure. First you map your geographic limits into matrix row and column limits and construct vectors of row and column indices to be read. Then you read the submatrix using readmtx, and construct the maplegend.
% we want to read latlim = [41 44]; lonlim = [-72 -69]; scalefactor = 1; % file format and contents nrows = 6293; ncols = 13461; lat1 = 71.36249997; lon1 = -179.13749993; yperrow = -dms2deg(mat2dms(0,0,30)); xpercol = dms2deg(mat2dms(0,0,30)); % map latlim and lonlim to row and column limits [rlim,clim] = yx2rc(latlim,lonlim,lat1,lon1,yperrow,xpercol) rlim = sort(rlim); readrows = rlim(1):scalefactor:rlim(2); readcols = clim(1):scalefactor:clim(2); % extract the map matrix map = flipud(readmtx('usa.bil',... nrows,ncols,'int8',readrows,readcols)); % construct maplegend cellsize = scalefactor*xpercol; maplegend = [1/cellsize, latlim(2)-cellsize/2, ... lonlim(1)-cellsize/2]; % display result usamap(latlim,lonlim,'lineonly') setm(gca,'mapproj','mercator'); tightmap meshm(map,maplegend)
Boston appears to have grown since the data in the VMAP0 was compiled. The VMAP0 also missed the heavily populated part of Cape Cod.
Some data is stored on regular grids in projected coordinates. In that case, you need to define a map projection that matches the one used for the data. Because the grid is in projected coordinates rather than latitude and longitude, the map is a general matrix map and needs a graticule. To construct one, build a matrix of row and column indices, convert them to projected coordinates, and carry out the inverse projection. The resulting latitude, longitude, and map matrices can be displayed in any projection. Edit the avhrrlambert interface function to see an example.
Another common type of geographic data file format is fixed length records with fields of information. Examples of such data include geographic names, hydrographic soundings, ship track lines, and astronomical catalogs of stars. Use the readfields function to read selected fields or records from such files.
As an example, retrieve the U. S. Geological Survey Geographic Names Information System (GNIS) file containing a concise list of names and locations within the United States. The file can be found at<http://mapping.usgs.gov/pub/gnis/US_concise.gz>.
The README file describes the file format and contents. To read this data, create a structure with this information.
filestruct(1).length = 51; filestruct(1).name = 'Feature Name'; filestruct(1).type = 'char'; filestruct(2).length = 10; filestruct(2).name = 'Feature Type'; filestruct(2).type = 'char'; filestruct(3).length = 32; filestruct(3).name = 'County'; filestruct(3).type = 'char'; filestruct(4).length = 17; filestruct(4).name = 'State'; filestruct(4).type = 'char'; filestruct(5).length = 16; filestruct(5).name = 'Geographic Coordinates'; filestruct(5).type = 'char'; filestruct(6).length = 6; filestruct(6).name = 'Feature Elevation'; filestruct(6).type = 'char'; filestruct(7).length = 1; filestruct(7).name = ''; %line ending filestruct(7).type = 'char';
You can identify records that match a search string using the grepfields function. With no output arguments, records matching the string are displayed onscreen. When an output argument is provided, the record numbers of matching records are returned.
Search for the locations of The MathWorks and Systems Planning and Analysis.
grepfields('us_concise','natick') Natick ppl Middlesex Massachusetts 421700N0712100W 180 indx = grepfields('us_concise','alexandria');
Use the readfields function to read the selected records into a structure.
s = readfields('us_concise',filestruct,indx) s = 1x12 struct array with fields: FeatureName FeatureType County State GeographicCoordinates FeatureElevation s(7) ans = FeatureName: 'Alexandria' FeatureType: 'ppl' County: 'Independent City' State: 'Virginia' GeographicCoordinates: '384817N0770250W' FeatureElevation: '32'
If you are having trouble getting the file format right, you can check the record length of a file with line ending characters. Remember that DOS-based computers use two line ending characters.
fid = fopen('us_concise','r'); str = fgets(fid) str = Aasu ppl Western American Samoa 141751S1704530W recordlength = length(str) recordlength = 133.00 fclose(fid);
If the file does not have line endings, you can use fread to read a few characters at a time.
| | Astronomical Data | Reference | |
