[GIS] Lookup Projection Name from EPSG with pyproj

projpyprojpython

Is it possible to programmatically lookup a projection's name (e.g. Albers Equal-Area) from its EPSG using pyproj?

I see pyproj.pj_list, but it does not contain a mapping to the EPSG.

Best Answer

List pj_list gives you at least projections the proj command understands. The EPSG context in the proj.4 biotop is a simple file based "DATABASE" that connects ONLY EPSG number entries (KEYS) with initialization params of the tool proj. You will find the data in the file 

 /usr/share/proj/epsg

if you use proj.4 and one of it's bindings in LINUX.

Each entry has an EPSG key and a set of proj params. It looks like this. A line beginning with # is a comment line and a line beginning with <NUMBER> is an EPSG data entry. Some example lines:

$ less /usr/share/proj/epsg

Output (..shorted)

...
# NAD83 / Texas Centric Albers Equal Area
<3083> +proj=aea +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=1500000 +y_0=6000000 +datum=NAD83 +units=m +no_defs  <>
# NAD83(HARN) / Texas Centric Lambert Conformal
<3084> +proj=lcc +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=1500000 +y_0=5000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs  <>
# NAD83(HARN) / Texas Centric Albers Equal Area
<3085> +proj=aea +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=1500000 +y_0=6000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs  <>
# NAD83 / Florida GDL Albers
<3086> +proj=aea +lat_1=24 +lat_2=31.5 +lat_0=24 +lon_0=-84 +x_0=400000 +y_0=0 +datum=NAD83 +units=m +no_defs  <>
...

You can use LINUX command line programs like grep to retrieve your info by screening the parameter +proj=aea (aea = Albers Equal Area) and the comments above the EPSG definition.

grep -B 1 'proj=aea' /usr/share/proj/epsg

Output (..shorted)

# NAD27 / Alaska Albers
<2964> +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +datum=NAD27 +units=us-ft +no_defs  <>
--
# NAD83 / BC Albers
<3005> +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs  <>
--
# NAD83 / Texas Centric Albers Equal Area
<3083> +proj=aea +lat_1=27.5 +lat_2=35 +lat_0=18 +lon_0=-100 +x_0=1500000 +y_0=6000000 +datum=NAD83 +units=m +no_defs  <>
--

And the parameter meaning of grep is:

grep -B 1 PATTERN FILE
      |       |     |
      |       |     +--- the file you want to screen ../usr/share/proj/epsg
      |       +--- the search pattern ..a line containing 'proj=aea'
      +--- show me one line before the pattern match happened ..usually the comment

You can write a small wrapper in python if you SHELL OUT the grep stuff in combination with the pj_list and build your own dictionary. Here a small "dirty" example (..for the shell out part only):

#!/usr/bin/python
import os

# Command
cmd  = 'grep -B 1 "proj=aea" /usr/share/proj/epsg'

# Shell out and handler
hnd  = os.popen(cmd)

# line to parse
ln   = "--"

# variable for comment (line before)
cmt  = ""

# variable for epsg definition
epsg = ""

# dict of found stuff
found = {}

# iterate lines in handler
while ln:

  # real on line
  ln = hnd.readline()

   # have a comment
  if ln.startswith('#'):
      cmt = ln

  # have a epsg definition 
  if ln.startswith('<'):

      # grab only the key 
      epsg = ln.split(' ')[0]

      # store in dict
      found[epsg]=cmt

# close handle after iteration 
hnd.close()

# show the results    
for epsg in sorted(found):
    print epsg, found[epsg]
Related Question