Every now and then I get the urge to define my own projection. Usually, I sit down for a while, hit my head on the wall, and the urge passes. For a few years I have worked with the Lake Erie and Allegheny Partnership for Biodiversity on various projects. Now we are getting deep into region-wide data collection, and so I decided to define an Albers Equal Area projection for the task. Specifically, I sought to refine an existing projection to match the bounds and center of the LEAP region. Yes, I know. I will be cursed one day for this decision, but it sure beats the alternatives at the moment.

Beyond defining the projection, I wanted it as an ESRI WKT and in Proj4 format. Here are the steps I took. I used spatialreference.org to help me do the format translation.

http://spatialreference.org needs the projection input in a form called Well Known Text (WKT)– specifically the Open Geospatial Consortium’s for of WKT. First, I uploaded the description in ESRI’s WKT: this doesn’t work. So here’s ESRI’s WKT for the LEAP boundary as I created a few days ago:

PROJCS["LEAP_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-80.75],PARAMETER["Standard_Parallel_1",39],PARAMETER["Standard_Parallel_2",43],PARAMETER["Latitude_Of_Origin",41],UNIT["Meter",1.0]]

I changed 4 things from the North America version of this:

"Central_Meridian"
"Standard_Parallel_1"
"Standard_Parallel_2"
"Latitude_Of_Origin"

OGC WKT has different names for these, specifically:

"longitude_of_center"
"Standard_Parallel_1"
"Standard_Parallel_2"
"latitude_of_center"

So, we can take USA Contiguous Albers Equal Area Conic from http://spatialreference.org/ref/esri/102003/

Go to it’s OGC WKT page:

http://spatialreference.org/ref/esri/102003/ogcwkt/

And thus extract the OGC WKT we want to modify:

PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",-96],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["latitude_of_center",37.5],UNIT["Meter",1],AUTHORITY["EPSG","102003"]]

Modifying those parameters thusly:

PARAMETER["longitude_of_center",-80.75],
PARAMETER["Standard_Parallel_1",39],
PARAMETER["Standard_Parallel_2",43],
PARAMETER["latitude_of_center",41]

Resulting in the following:

PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["longitude_of_center",-80.75],PARAMETER["Standard_Parallel_1",39],PARAMETER["Standard_Parallel_2",43],PARAMETER["latitude_of_center",41],UNIT["Meter",1]]

Which I have uploaded to a custom LEAP projection.

Ok. Final step. gdalwarp was the original tool I wanted to use. It requires that we define our projection in a format called Proj4 or using an EPSG code. Since we invented this projection, it’s has no EPSG code. Now that we have a definition for it loaded into http://spatialreference.org, we can use the web site to give us the proj4 definition:

http://spatialreference.org/ref/sr-org/7915/proj4/

This site now gives us the following:

+proj=aea +lat_1=39 +lat_2=43 +lat_0=41 +lon_0=-80.75 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs

The sample reprojection code from the gdalwarp website is as follows:

gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif

To use it for our own data, we’ll do something like this:

gdalwarp -t_srs '+proj=aea +lat_1=39 +lat_2=43 +lat_0=41 +lon_0=-80.75 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ' input.tif output.tif

Boom! Data reprojected. And so another projection is born. Oops.

postscript– Yes I know about epsg.io— as far as I can tell, it does not yet have the ability to receive uploads.

Edit from Howard Butler:

@smathermather if your GDAL has curl support, you can give http://t.co/C5evFSK352 URLs to custom SRSs and it'll just use 'em.

— Howard Butler (@howardbutler) March 18, 2014

ala:

gdalwarp -t_srs "http://spatialreference.org/ref/sr-org/7915/proj4/" input.tif output.tif

Brilliant!