Learning the QGIS Python API - Building Mapping Applications with QGIS - page 77

Analyzing raster data
We're going to start by writing a program to load in some raster-format data and analyze its
contents. To make this more interesting, we'll use a Digital Elevation Model ( DEM ) file,
which is a raster format data file that contains elevation data.
The Global Land One-Kilometer Base Elevation Project ( GLOBE ) provides free DEM
data for the world, where each pixel represents one square kilometer of the Earth's surface.
GLOBE data can be downloaded from http://www.ngdc.noaa.gov/mgg/topo/gltiles.html .
Download the E tile, which includes the western half of the USA. The resulting file, which
is named e10g , contains the height information you need. You'll also need to download the
e10g.hdr header file so that QGIS can read the file—you can download this from ht-
tp://www.ngdc.noaa.gov/mgg/topo/elev/esri/hdr . Once you've downloaded these two files,
put them together into a convenient directory.
You can now load the DEM data into QGIS using the following code:
registry = QgsProviderRegistry.instance()
provider = registry.provider("gdal", "/path/to/e10g")
Unfortunately, there is a slight complexity here. Since QGIS doesn't know which coordin-
ate reference system is used for the data, it displays a dialog box that asks you to choose
the CRS. Since the GLOBE DEM data is in the WGS84 CRS, which QGIS uses by default,
this dialog box is redundant. To disable it, we need to add the following to the top of our
program:
from PyQt4.QtCore import QSettings
QSettings().setValue("/Projections/defaultBehaviour",
"useGlobal")
Now that we've loaded our raster DEM data into QGIS, we can analyze it. While there are
lots of things we can do with DEM data, let's calculate how often each unique elevation
value occurs within the data.
Note
Notice that we're loading the DEM data directly using QgsRasterDataProvider . We
don't want to display this information on a map, so we don't want (or need) to load it into
QgsRasterLayer .