Mapping with Python
This page outlines the pythonic solution I developed when I had to make several maps of the same spatial domain but with different data. I typically use QGIS and GRASS GIS but the process to output a final map takes 10 to 15 minutes and when producing 50 maps, that is too long. With what we cover here, the same task can be completed in moments.
- Environment
The first thing to address is the environment. I like to use Anaconda which makes this task quite straightforward. To create a new environment, navigate to where you want your project to “live”. Then open a command prompt (I prefer cmder) and enter…
conda create — name YourEnvironmentName python=3.6
This will usually take a little while to complete and it may ask if you would like to proceed. Enter ‘y’ to indicate you would like to proceed and once it has finished, you can activate the environment.
activate YourEnvironmentName
This should activate your environment. You can tell because the cursor should have YourEnvironmentName in front of it. Note that sometimes you have to (or you can) use
conda activate YourEnvironmentName
You will have to install some packages. They can be installed using pip or conda. Conda tends to be more effective.
conda install basemap
If you are curious about where your newly created virtual environment “lives”, the files for all of your virtual environments should be at C:\Users\~user~\Anaconda3\envs\ on a windows system. We should now be set up for mapping! If you have any trouble with installation, feel free to reach out and I will try to help. I have done this with both a Python 2.7 and 3.6 environment.
2. Imports
Now let’s get into it! I am a scripter so I will just go through my script piece by piece. I like to start with all of my imports so that I have it when I need it. I will continue to display text intended to be code with the command line style formatting so as to be clear about what is code and what is explanation.
from mpl_toolkits.basemap import Basemap, pyproj
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from pylab import *
import os
import glob
3. Functions
Next I like to import my functions. Again, this is advantageous to do early in the script because it can be called anytime afterward. Some of these are pretty straightforward. Some are just part of how I code. Others are very specific to the task at hand, mapping. I’ll point out and explain what they are each for.
# Define last position finder
def find_last(s,t):
last_pos = -1
while True:
pos = s.find(t, last_pos +1)
if pos == -1:
return last_pos
last_pos = pos
This is a pretty simple function that I mainly use to navigate the path variable. I will usually call the relative path to my script and then, because I use a convention of keeping my Python files in a \Python folder for projects, I need to get out of that folder in my path to access other files. Hopefully that makes more sense for you shortly.
4. File Structure
I like to write scripts that run in any place and create their own folders if needed. For that, we have to do some work with the path and make sure that your file structure is correct. The image below shows how you should set up your file structure.
The .py file should go in the Python folder. Download this static file of the New Mexico border and move it into the Static Files folder. This is the file that will be the same on each map. Then download these dynamic data files ( NM Counties, NM State Parks, BLM Wilderness Study Areas, Land Grants, American Indian Lands ) and move them to the Data folder.
With that said, let’s jump right into the code.
# Define path
path = os.path.abspath(os.path.dirname(__file__))
# Shorten path to one folder up
path = path[:find_last(path,os.sep)]
# Unzip the zipped static shp files
for file in glob.glob(path+os.sep+’Static Files’+os.sep+’*.zip’):
zip_ref = zipfile.ZipFile(file, ‘r’)
zip_ref.extractall(path+os.sep+’Static Files’+os.sep)
zip_ref.close()
# Get the name of the static file
staticFile = glob.glob(path+os.sep+’Static Files’+os.sep+’*.shp’)
staticFile = staticFile[0]
# Unzip the zipped dynamic shp files
for file in glob.glob(path+os.sep+’Data’+os.sep+’*.zip’):
zip_ref = zipfile.ZipFile(file, ‘r’)
zip_ref.extractall(path+os.sep+’Data’+os.sep)
zip_ref.close()
# Create array to store names
Files = []
# Get list of files to map
for file in glob.glob(path+os.sep+’Data’+os.sep+’*.shp’):
Files.append(file)
Files.append(file)
Here we get the .py file’s path and use it to find the other files to map. The .shp files that we downloaded are zipped and need to be unpacked to be accessed by the Basemap module so this bit of code does that. It then also builds a list of the files that will be mapped on each seperate figure. Go ahead and try for yourself. You should have something similar to this.
>>>path
‘C:\\Tools\\Basemap Mapper’
>>>Files
[‘C:\\Tools\\Basemap Mapper’\\*.shp’,
‘C:\\Tools\\Basemap Mapper’\\*.shp’,
‘C:\\Tools\\Basemap Mapper’\\*.shp’]
It’s unlikely that your file name will print out so cleanly but you should have a list of however many files you have in the Data folder, 5 files if you downloaded all of those referenced on this page. You can check the length to be sure. You can also check the staticFile to be sure that it is a single string representing the path to the file. That should have us ready to start the actual mapping now!
5. Mapping
Here is the whole code. I’ll break it down a little below.
# Loop through files
# — — — — — — — — — — — — — — — — — — — — — — — — — — —
for file in Files:
# Define the bounds of the map area
LongLeft = -110
LongRight = -102.5
LatBottom = 31
LatTop = 38
midLat = (LatBottom+LatTop)/2.0
midLong = (LongLeft+LongRight)/2.0
# Open a figure for plotting
plt.figure(figsize=(8, 9))
# Create the map object
map = Basemap(
llcrnrlon=LongLeft,
llcrnrlat=LatBottom,
urcrnrlon=LongRight,
urcrnrlat=LatTop,
projection=’tmerc’, lat_0 = midLat, lon_0 = midLong, epsg=3857)
# Get imagery
map.arcgisimage(service=’ESRI_Imagery_World_2D’, xpixels = 1500, verbose= True)
# Read in the static shp files
map.readshapefile(staticFile[:-4], staticFile[find_last(staticFile,os.sep)+1:-4], drawbounds = True, color=’k’, linewidth=3)
# Get the dynamic file name for naming and handling
fileName = file[find_last(file,os.sep)+1:-4]
# Open data file
dynamicFile = map.readshapefile(file[:-4], fileName, drawbounds = True, color=’w’, linewidth=1.5)
# Make a title
plt.title(fileName, fontname=’Times New Roman’, y=1.08, fontsize=14, color= ‘k’, fontweight=’bold’)
# Check to make sure directory exists for saving figure
if not os.path.exists(path+os.sep+’Output’+os.sep):
os.makedirs(path+os.sep+’Output’+os.sep)
# Save figure
plt.savefig(path+os.sep+’Output’+os.sep+fileName+’.png’,dpi=500)
plt.close()
First we open a loop to move through the dynamic files that we want to plot. Then we define the domain that we would like to display. We also define the middle of the domain to be passed to the Basemap module object when we call that. Next we open a figure to plot in.
Then we create the map object by calling the Basemap function* from the basemap module. This takes in parameters for the lower left corner (llcrnrlon & llcrnrlat) and the upper right corner (urcrnrlon & urcrnrlat). It takes these inputs in latitude and longitude. It then allows the user to define the projection. I used the ‘tmerc’ or the Transverse Mercator Projection in which the globe is first rotated so the central meridian becomes the “equator”, and then the normal mercator projection is applied. The middle of the domain is supplied and then the EPSG code is defined. We use 3857 here. For more information on this EPSG code, click here.
Next, we get the aerial imagery. This loads a basemap for whatever bounds you have defined. After this, we open the static .shp file that you want displayed on every map. We have the border of New Mexico for this since that is the theme of the maps. This file will be the same each time. The next step is to get the dynamic files that we are looping through.
First, we get the name from the ‘file’ item that we are looping through and use it to name the map object that we are adding to the figure. It is then used to add a title to the figure. Finally, we save the figure after checking that the ‘Output’ folder exists to save figures into. We then close the figure and start the loop over.
This process produces the following figures. You can use the process to map whatever .shp files you need or want to.
How can I help you?
I have a master’s of science in hydrogeology. In addition to technical environmental skills, I have a passion for computer programming and solving challenging computer based problems. I have worked with huge data sets, worked with and developed databases with myriad data types, written basic software to work around proprietary software limitations, automated rote tasks, automated complex tasks…
I would love to talk to you about my products and services and how we might work together.
If you would like to get in touch with me to discuss a potential project or to engage my services, please take action and email or call me now
