Python | Geospatiality

This blog has moved to a new home at Final Draft Mapping.

Generate a Projection File (.prj) using Python

A .prj file contains the coordinate system information for a dataset and is required for ‘on the fly’ projection by GIS software. The file itself is a text file containing information in Well Known Text (WKT) format. The code snippet here shows you how to generate a .prj file, using Python, for your data if the .prj file is missing. See the post on CSV to Shapefile with pyshp for an example of using the function in a workflow. You will need to know what coordinate system the data is in and the EPSG code. For example, to find the EPSG code for WGS84 use a search engine and search for EPSG WGS84. The first site returned is usually http://spatialreference.org/. The EPSG code for WGS84 is 4326.

# function to generate .prj file information using spatialreference.org
def getWKT_PRJ (epsg_code):
     import urllib
     # access projection information
     wkt = urllib.urlopen("http://spatialreference.org/ref/epsg/{0}/prettywkt/".format(epsg_code))
     # remove spaces between charachters
     remove_spaces = wkt.read().replace(" ","")
     # place all the text on one line
     output = remove_spaces.replace("\n", "")
     return output

# create the .prj file
prj = open("filename.prj", "w")
# call the function and supply the epsg code
epsg = getWKT_PRJ("4326")
prj.write(epsg)
prj.close()

The filename should match the name of your data files. For examples, the .prj file for cities.shp or cities.tab should be cities.prj.

If you open the .prj file in a text editor you will see the text below.

GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

GEOGCS at the beginning of the WKT indicates that EPSG: 4326 is a geographic coordinate reference system (it would be PROJCS if it was a projected coordinate system). The information within the square brackets after DATUM provide the information for the parameters of the datum. We see the name of the datum, ‘WGS_1984’, and the SPHEROID information with a semimajor axis of 6,378,137 m with an inverse-flattening ratio of 298.257223563. The PRIMEM data tells us that it uses Greenwich as the prime meridian (where longitude is zero). UNIT specifies the measurement unit of the coordinate system, here it is ‘degree’, and the 0.0174532925199433 value is required to convert from radians into the stated units.

Open your data in a GIS and make sure that it is displaying correctly in relation to other features.

import math def dd2dms(longitude, latitude): # math.modf() splits whole number and decimal into tuple # eg 53.3478 becomes (0.3478, 53) split_degx = math.modf(longitude) # the whole number [index 1] is the degrees degrees_x = int(split_degx[1]) # multiply the decimal part by 60: 0.3478 * 60 = 20.868 # split the whole number part of the total as the minutes: 20 # abs() absoulte value - no negative minutes_x = abs(int(math.modf(split_degx[0] * 60)[1])) # multiply the decimal part of the split above by 60 to get the seconds # 0.868 x 60 = 52.08, round excess decimal places to 2 places # abs() absoulte value - no negative seconds_x = abs(round(math.modf(split_degx[0] * 60)[0] * 60,2)) # repeat for latitude split_degy = math.modf(latitude) degrees_y = int(split_degy[1]) minutes_y = abs(int(math.modf(split_degy[0] * 60)[1])) seconds_y = abs(round(math.modf(split_degy[0] * 60)[0] * 60,2)) # account for E/W & N/S if degrees_x < 0: EorW = "W" else: EorW = "E" if degrees_y < 0: NorS = "S" else: NorS = "N" # abs() remove negative from degrees, was only needed for if-else above print "\t" + str(abs(degrees_x)) + u"\u00b0 " + str(minutes_x) + "' " + str(seconds_x) + "\" " + EorW print "\t" + str(abs(degrees_y)) + u"\u00b0 " + str(minutes_y) + "' " + str(seconds_y) + "\" " + NorS # some coords of cities coords = [["Dublin", -6.2597, 53.3478],["Paris", 2.3508, 48.8567],["Sydney", 151.2094, -33.8650]] # test dd2dms() for city,x,y in coords: print city + ":" dd2dms(x, y)

Geospatiality

make geospatial your speciality

Tag Archives: Python

Generate a Projection File (.prj) using Python

Decimal-Degrees (DD) to Degrees-Minutes-Seconds (DMS) with Python