Making static maps and processing geodata with GIS

Introducing QGIS

GQIS is the leading free, open source Geographic Information Systems (GIS) application. It is capable of sophisticated geodata processing and analysis, and also can be used to design publication-quality data-driven maps.

The possibilities are almost endless, but you don’t need to be a GIS expert to put it to effective use for both displaying and processing geographic data, whether for online or print maps.

Launch QGIS, and you should see a screen like this:

The data we will use today

Download the data from this session from here, unzip the folder and place it on your desktop. It contains the following folders and files:

Map seismic risk and earthquakes

Make a choropleth map showing seismic risk for the continental United States

As we discussed in week 9, choropleth maps fill areas with color, according to values in the data.

We will first import the shapefile in the seismic folder.

Select Layer>Add Layer>Add Vector Layer... or click on this icon:

At the dialog box, click Browse and navigate to the file seismic. It is important that you select the file with the .shp extension. Then click open, and the following map should appear, filled with a random color:

The Layers panel ay bottom-left should now show the seismic layer. You can turn off the visibility of any layer by unchecking its box in the Layers panel. This can be useful to see the status of layers that would otherwise be obscured.

These controls allow you to pan and zoom the display:

You can focus the display on the full extent of any layer by right-clicking it in the Layers panel and selecting Zoom to layer.

Notice EPSG:4326 at the bottom right. This defines the map projection and datum for the layer.

Right-click on seismic in the Layers panel at left and select Properties>General. You should see the following under Coordinate reference system:

EPSG:4326 and WGS84 mean that this shapefile is in an equirectangular projection, and the WGS84 datum. This is also the default for QGIS if no projection is specified.

We will select another projection for our map later. Click Cancel or OK to close Properties for this layer.

Now we need to color the areas for the sesimic layer by values in the data. Right-click on the layer in the Layers panel and select Open Attribute Table, which corresponds to the .dbf of the shapefile:

There is one variable, ValueRange, which gives ranges for the percentage chance of experiencing a damaging quake in 2017.

Close the attribute table and open Properties>Style for the seismic layer. Select Categorized from the dropdown menu at top, which is the option to color data according to values of a categorical variable.

(Note, although these values are percentages, they have already been binned into categories in this data. If we were coloring by a continuous variable that had not already been binned, we would select Graduated and select values for the bins. See last year’s class for an example.)

Click the Classify button and the values in the data will be assigned random colors:

Notice that the range 10-12 is out of order, and we also have a blank range at the bottom. First select this blank and click the red minus symbol to remove it. Now click and drag the range 10-12 to the correct position.

Now click on the Color ramp dropdown menu and select Reds:

It is also possible to import a wide variety of color palettes, including those from ColorBrewer, using the Color ramp dropdown menu.

Click OK, and the map should look like this:

Next we’ll remove the black outlines around the polygons. Open Properties>Style again, and double-click on the colored square under Symbol. Select Simple fill in the dialog box:

Scroll down, and change Outline style to No Pen:

Click on OK at each of the dialog boxes and the map should now look like this:

It is also possible to edit the colors for each value individually.

If you are likely to want to style data with the same column headings in the same way in future, it is a good idea to click the Style button at bottom left and select Save Style>QGIS Layer Style File. This saves the style in QML, which is a variant of XML. When loaded using the Style>Load Style, it will automatically apply the saved styling to future maps.

Now is a good time to give the project a projection: We will use a variant of the Albers Equal Area Conic projection, optimized for maps of North America.

Select Project>Project Properties>CRS (for Coordinate Reference System) from the top menu, and check Enable 'on the fly' CRS transformation. This will convert any subsequent layers we import into the Albers projection, also.

Type Albersinto the Filter box and select North_America_Albers_Equal_Area_Conic, which has the code EPSG:102008:

Click OK and the map should reproject. Notice how EPSG:102008 now appears at bottom right:

Save the project, by selecting Project>Save from the top menu.

Add a base layer showing countries

Add the ne_50m_admin_0_countries_lakes layer, as before. It will appear above our choropleth map:

So click and drag the seismic layer over the base layer in the Layers Panel. Also right-click on the layer and Zoom to Layer so that the view is optimized for this projection:

Now open up the Properties for the ne_50m_admin_0_countries_lakes layer, select Simple fill, and click on the Fill color to change it to #cccccc for a light gray. Also change the Outline to #ffffff for white:

Click OK and the map should look like this:

Now select Project>Project Properties...>General, click on the Background color and set it to #686868 for a dark gray:

The map should now look like this:

Add the earthquakes layer

To import a CSV or other delimited text file with points described by latitude and longitude coordinates, select Layer>AddLayer>Add Delimited Text Layer from the top menu or click this icon:

Browse to the file quakes.csv, and ensure that the dialog box is filled in like this:

If your text file is not a CSV you may have to select the correct delimiter, and if your latitude and longitude fields have other names, you may have to select the X field (longitude) and Y field (latitude) manually.

When you click OK you will be asked to select a projection, or CRS, for the data. You may be tempted to select the same Albers projection we have set for the project, but this will cause an error. QGIS will handle the conversion to that projection: Because the data in the CSV file is not yet projected, instead select a datum with the default equirectangular projection WGS 84 EPSG:4326 (NAD83 EPSG:4269, which is an equirectangular projection in the 1983 North American Datum, would also work.)

Click OK and a large number of points will be added to the map:

Now we will style these points, giving them a single color but sizing them according to the amount of shaking that quake each caused. Open the attribute table for the quakes layer and notice that it contains a variable called mag, for earthquake magnitude. This is a logarithmic scale, so that a magnitude difference of 1 corresponds to a 10-fold difference in earth movement, as recorded on a seismograph, as this calculator shows:

This means that the equation to convert magnitude into the amplitude of the shaking, as recorded on a seismograph, is amplitude = 10^magnitude, or ten raised to the power of magnitude.

So if we want to scale the points according to the actual size of the quake, we need to create a new field in the data for amplitude. To do this, click on the abacus symbol in the top toolbar:

Fill in the dialog box as follows:

The formula is: amplitude = 10^mag/100000. The final division by 100,000 helps keep the values in a range that is easy to read and understand, which will help when we set the size of the points. The Outout preview at the bottom shows the resulting value for the first row in the attribute table. Click OK to create the new amplitude variable.

Select Properties>Style for the quakes layer, set the Color to #ffffff for white and the transparency to 50%:

Now click on this icon for Size, and select Size Assistant...:

Fill in the dialog box as follows:

Using the Scale method of Surface will scale the circles by area, according to values in the data; Flannery is a method of scaling that is designed to compensate for the errors people make when estimating the area of circles. You may wish to experiment with this, but I would strongly advise against using Radius, for the reasons we discussed in week 2.

By default, QGIS will select the range of values in the data. Make sure to set both Size from and Value from to zero to ensure that the area of the points scales accurately according to values in the data.

Click OK until you return to the map, which should now look like this:

Export the finished map in raster image or vector graphic formats

To export the map as a JPG, PNG, or other raster image format, select Project>Save as Image... from the top menu.

To export as a vector graphic, for further editing in Adobe Illustrator or another vector graphic editing program, select Project>New Print Composer, give the composer an appropriate name and click OK. Set Page size as required, then click the Add new map icon:

Draw a rectangle over the page area, and the map should appear. At this point, you may need to alter the zoom level and position of your map in the main QGIS display to get a pleasing display in the print composer.

Your print composer should now look something like this:

The print composer also allows you to add and edit a legend, by clicking this icon:

You can also add text and other items. In practice, however, you may find it easier to add legends and labels in Adbobe Illustrator.

The map can be exported in SVG and PDF vector formats by clicking these export icons:

Finally, save the QGIS project.

Join external data to a shapefile

Open a new project by selecting Project>New from the top menu.

In this new project, import the ne_50m_admin_0_countries_lakes shapefile. Right-click on it in the Layers panel and Save As ... an ESRI shapefile, keeping the default WGS 84 CRS. Browse to the gdp_pc folder, call the new file gdp_pc and check the option to Add saved file to map.

Now remove the original shapefile, by right-clicking on that layer and selecting Remove.

Open the attribute table of the new shapefile and notice that it contains a field called iso_a3, which is a three-letter code for each country, assigned by the International Organization for Standardization.

Now use Add Vector Layer to import the file gdp_pc.csv. (Note that when joining external data in a CSV file to a shapefile, you do not import the file as a delimited text file, as we did previously to display data on a map.)

After import, this file will appear as an isolated table in the Layers panel.

Select it and open the attribute table to view the data, which contains country names, the three-letter country codes, and data on GDP per capita in 2015:

Notice that some cells contain the value -99, which is here being used to designate null values, where there is no data.

The same subfolder with the file also contains the file gdp_pc.csvt. This contains information about the type of data in each field of the CSV file, in this case:

"String","String","Real"

When we make the join, this will tell QGIS what sort of data is in each field of the CSV file. String indicates a string of text, Real indicates numbers that can include decimals, while Integer indicates whole numbers. Without this information, QGIS will treat all of the fields in the file as text.

Close the attribute table once more, right-click on the gdp_pc shapefile and select Properties>Joins. Click the green plus sign and fill in the dialog box as follows to join the CSV file to the shapefile by the iso_a3 three-letter country codes:

Click OK at each dialog box to completethe join, then open shapefile’s attribute table once more to confirm that the data from the CSV file has appeared.

Save the joined data in another geodata format

Right-click on the joined shapefile, select Save As ... and notice the Format options include ESRI shapefile, GeoJSON and KML. You are also able to choose a projection (CRS) for the new file and restrict its extent by latitude and longitude coordinates. Click the Browse button to select the location to save the data.

Save this file in the gpd_pc folder as GeoJSON with an appropriate name, keeping the default WGS 84 CRS.

Use QGIS’s vector geoprocessing tools

Start a new project, and import both the ne_50m_admin_0_countries_lakes shapefile, and the seismic_raw shapefile. Note that it extends beyond the boundaries and coastline of the U.S.:

Turn off the visibility of the seismic_raw layer, turn on the visibility of the ne_50m_admin_0_countries_lake layer, select it in the Layers Panel, and click the Select feature(s) button:

Now click anywhere within the borders of the United States to select the nation, which should now appear yellow:

Select Vector>Geoprocessing Tools>Intersection from the top menu and fill in the dialog box as follows:

Click Run and a new layer will be created, clipped to the borders and coastline of the United States:

Its attribute table will contain data from both layers used for the intersection.

You can now right-click on the new Intersection layer and Save As ... as a shapefile, GeoJSON, or another geodata format.

Look at the other options in the Vector>Geoprocessing Tools menu. Their icons give a good idea of what they do. (Clip is similar to Intersection, except that it does not include data from both layers in the new attribute table.)

Now open a new project, and import the shapefile sf_test_addresses. I made this shapefile from the addresses we geocoded in week 9. I saved it in a Google Mercator projection, also known as EPSG:900913, used for Google and other online maps.

This is important, because we are going to create a “buffer” defining areas within 1,000 feet of the nearest point. For this we need a projection with units defined in distance, rather than degrees, which is the unit for the WGS 84 datum.

To confirm the units for the Google Mercator projection, click on the globe symbol at bottom-right:

You should see that the CRS/projection information that appears near the bottom of the dialog box contains +units=m, which tells us that distances in this projection are measure in meters:

Creating a buffer shapefile is a task you might perform if, for example, working out which areas are off-limits for sex offenders under residency restrictions.

Select Vector>Geoprocessing Tools>Fixed distance buffer and fill in the dialog box as follows:

Selecting 99 under Segments to approximate ensures that the resulting shapes are as smooth circles. Buffer distance is set to 304.8 because the projection’s units are meters; this value gives the 1,000 feet we require. Checking Dissolve result merges overlapping buffers into the same polygon.

Click Run, and this should be the result:

Again, you can now right-click on the new Buffer layer and Save As ... as a shapefile, GeoJSON, or another geodata format.

Assignment

File an update on progress with your final project. Pay particular attention to identifying problems/obstacles.

You can file the update as a markdown file in GitHub, or by email. Either way, make sure that your GitHub repo is up to date with data, charts etc, so that I can sync with my version of your project and review your progress. The update should be sufficiently detailed that I can look at your repo and understand what you are doing.

If you are having trouble with GitHub, get in contact for a refresher.

Further reading

QGIS Training Manual

QGIS User Guide

A Gentle Introduction to GIS