Python: Accessing GIS Data With URL and FTP Modules

Hello Readers,


One of the most utilized methods in obtaining GIS data in Python is through URL and FTP connections. This post will demonstrate reading file data from the U.S. Geological Survey (USGS) and the National Oceanic and Atmospheric Administration (NOAA).

URL

The Uniform Resource Locator enables web browsers and programs to locate resources. We will use the urllib module to access the file, "hancock.zip", from the Google code website. After opening the Python Interpreter, or writing a script, import the urllib module, and assign the URL and file name to variables. After retrieving the file, Python will print out confirmation:

Retrieving hancock.zip

Now we turn to the USGS, and we want to read in US earthquake data from the past hour into Python. We use the urlopen() method to establish the connection. Using readline() to parse and print each line in the CSV file, we can see the first line includes the column headers, and the second line describes an earthquake near The Geysers in California. It appears to be small in magnitude at 0.7.

Reading Earthquake Data

Of course we can loop through the rest of the lines with a for loop:

More Earthquake Data

As we can see, the recently reported earthquakes all happened in the Western U.S.- Alaska, California, and Oregon.

FTP

The File Transfer Protocol allows file transfers between hosts over a TCP network. We can access NOAA's ftp server using the ftplib module to retrieve a specific file about tidal data. After specifying the server, we can connect to it, which requires no login information (anonymous). We change the directory and write the file in binary ("wb").

We look for a specific line concerning the location in latitude and longitude of this particular tidal buoy in the for loop. (It is located off the coast of Alaska- 50N, 171.8W)

Retrieving Tides Data via FTP

Likewise, we can approach this retrieval process using the urllib module as well. Since we can locate the tidal text file on the NOAA ftp server, we simply use "ftp://" as the URL protocol.

Retrieving Tides Data via URL

After running the same loop, we arrive with the same latitude and longitude.

And there we have it, folks! Two ways of retrieving and reading data from online resources through Python- yes, pure Python. Here we retrieved shapefiles from hancock.zip, earthquake data from USGS, and tidal buoy data from NOAA. Stay tuned for more posts!


Thanks for reading,

Wayne
@beyondvalence

Code:

# Python GIS URL FTP


# URL
import urllib

# access hancock shapefiles
url = "https://geospatialpython.googlecode.com/files/hancock.zip"
fileName = "hancock.zip"
urllib.urlretrieve(url, fileName)
print("\n Saved Hancock file \n")
# hancock.zip is saved to current working directory

# access earthquake data from USGS
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_hour.csv"
earthquakes = urllib.urlopen(url)
print("reading earthquake lines \n")
print(earthquakes.readline() + "\n")
print(earthquakes.readline() + "\n")

# iterate loop
print("looping earthquake data \n")
for record in earthquakes: print record

# FTP
import ftplib

server = "ftp.ngdc.noaa.gov"
dir = "hazards/DART/20070815_peru"
fileName = "21415_from_20070727_08_55_15_tides.txt"
ftp = ftplib.FTP(server)
ftp.login()
ftp.cwd(dir)
out = open(fileName, "wb")
ftp.retrbinary("RETR " + fileName, out.write)
out.close()
dart = open(fileName)
print("dart \n")
for line in dart:
    if "LAT, " in line:
        print line
        break

print ("dart URL \n")
# or use URL        
dart = urllib.urlopen("ftp://" + server + "/" + dir + "/" + fileName)
for line in dart:
    if "LAT, " in line:
        print line
        break
        
#
print("done")

Python: A Simple GIS Model of Colorado

Hello Readers,


Since I have not posted with Python in awhile, this is a great opportunity to explore geographic information systems (GIS) with Python! We shall start mapping with a simple graphic of the Colorado state with a few cities using turtle graphics.

Though we do not use fancy graphics here, the resulting map below introduces the GIS capability of pure Python. For further GIS work with a variety of data formats such as raster and vector data, we will use the geospatial data abstraction library (GDAL). However, with the simple turtle graphics, we can create the map showing the 3 cities in Colorado and their populations:

In addition to the city names and their populations, we also print the largest city and the western most city by longitude from querying the population and coordinates for the maximum and minimum, respectively. The largest city has the largest population, and the western most city has the lowest longitude, being negative in the western hemisphere.

The code is displayed below as a reference:

 # GIS in Python Part 1
 
import turtle as t
 
## Section 1: Data Model
NAME = 0
POINTS = 1
POP = 2
state = ["COLORADO", [[-109, 37],[-109, 41],[-102,41],[-102,37]], 5187582]
cities = []
cities.append(["DENVER", [-104.98, 39.74], 634265])
cities.append(["BOULDER", [-105.27, 40.02], 98889])
cities.append(["DURANGO", [-107.88, 37.28], 17069])
map_width = 400
map_height = 300
 
 # bounding box for COLORADO
minx = 180
maxx = -180
miny = 90
maxy = -90
for x,y in state[POINTS]:
    if x < minx: minx = x
    elif x > maxx: maxx = x
    if y < miny: miny = y
    elif y > maxy: maxy = y
dis_x = maxx - minx
dis_y = maxy - miny
x_ratio = map_width / dis_x
y_ratio = map_height / dis_y
 
def convert(point):
    lon = point[0]
    lat = point[1]
    x = map_width - ((maxx-lon) * x_ratio)
    y = map_height - ((maxy-lat) * y_ratio)
    # python turtle graphics start in middel of screen
    # must offset points to center
    x = x - (map_width/2)
    y = y - (map_height/2)
    return [x,y]
 
## Section 2: Map Renderer
# COLORADO
t.up()
first_pixel = None
for point in state[POINTS]:
    pixel = convert(point)
    if not first_pixel:
        first_pixel = pixel
    t.goto(pixel)
    t.down()
t.goto(first_pixel)
t.up()
t.goto([0,0])
t.write(state[NAME], align="center", font=("Arial", 16, "bold"))

# CITIES
for city in cities:
    pixel = convert(city[POINTS])
    t.up()
    t.goto(pixel)
    # place a point for the city
    t.dot(10)
    # label the city
    t.write(city[NAME] + ", Pop.: " + str(city[POP]), align="left")
    t.up

# attribute query, use keyword lambda
biggest_city = max(cities, key=lambda city:city[POP])
t.goto(0,-200)
t.write("The biggest city is: " + biggest_city[NAME])
western_city = min(cities, key=lambda city:city[POINTS])
t.goto(0,-220)
t.write("The western-most city is: " + western_city[NAME])
t.goto(0, 200)
t.write("MAP OF COLORADO", align="center", font=("Arial", 20, "bold"))
t.pen(shown=False)
t.done()

There will be more python posts with GDAL and colorful visualizations! Stay tuned, and in the meantime, take a look at 'greatest ever' infographic :



Thanks for reading,

Wayne
@beyondvalence