[Gdal-dev] fun with projections -- GCS to AEA

Fri Jun 16 11:22:15 EDT 2006

Hello all and happy Friday.  I have been working on something that isn't
coming out quite right so I though I would write it up for my own sanity
and share it in case someone was inclined to look it over.

I have a .shp of US counties in an unknown GCS.  QGIS assumes a WGS84,
good enough for me -- EPSG:4326 it is.  I need 50-mile buffers for each
county and my research tells me that whether I use GRASS or PostGIS I
need by base data to be projected into some cartesian units because the
buffer functions use map units.  Looks like a job for ogr2ogr.

But what projection?  Further digging indicated that an Albers
Equal-Area Conic is a good choice when working with the full extent of
the states.  I found some clues at (1), (2), and (3).  Neither OGR nor
QGIS seem to have any EPSG:9822 so I take this to be a general
transformation without the specific parameters.

So I wrote:

> g.proj -w proj4='+proj=aea +units=mi +datum=WGS84 +ellps=WGS84 \
    +lat_1=29.5 +lat_2=45.5 +lon_0=-96' > aeaUS.wkt

# aeaUS.wkt attached

> ogr2ogr -s_srs EPSG:4326 -t_srs aeaUS.wkt temp usagen counties
> ogrinfo -so temp counties > counties_aea.ogrinfo

# counties_aea.ogrinfo attached

What concerns me are the extents reported by ogrinfo:
Extent: (-6293474.061879, 2670863.569639) - (2256278.441196,
8557825.854853)

8.5 million miles wide?  I don't think so.

When I set a custom projection in QGIS using the same proj4 string and
add the original unprojected layer I see this metadata for the layer:

General:
Geometry type of the features in this layer : Polygon
The number of features in this layer : 3140
Extents:
In layer spatial reference system units : xMin,yMin -178.215,18.9248 :
xMax,yMax -66.9698,71.4066
In project spatial reference system units : xMin,yMin -4903.87,1195.6 :
xMax,yMax 1937.56,5684.6
Layer Spatial Reference System:
+proj=longlat +ellps=WGS84 +no_defs 
Project (Output) Spatial Reference System:
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=0 +lon_0=-96 +x_0=0 +y_0=0
+ellps=WGS84 +datum=WGS84 +units=mi

~ 6800 miles wide -- that's more like it.  Cook County, IL (Chicago) is
about 50 miles north to south.  It's about 160 miles from Chicago to
Indianapolis as the crow flies.  I can live with that, but I need a
projected layer on disk, not just QGIS doing it on the fly.

When I open the ogr-generated .shp in a new QGIS project I get 

General:
Geometry type of the features in this layer : Polygon
The number of features in this layer : 3140
Extents:
In layer spatial reference system units : xMin,yMin -6.29347e
+06,2.67086e+06 : xMax,yMax 2.25628e+06,8.55783e+06
In project spatial reference system units : xMin,yMin -6.29347e
+06,2.67086e+06 : xMax,yMax 2.25628e+06,8.55783e+06
Layer Spatial Reference System:
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=0 +lon_0=-96 +x_0=0 +y_0=0
+ellps=WGS84 +datum=WGS84 +units=mi
Project (Output) Spatial Reference System:
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=0 +lon_0=-96 +x_0=0 +y_0=0
+ellps=WGS84 +datum=WGS84 +units=mi

Cook County is now over 77,000 miles long!

Q: Is this a good choice of projections?

Q: What is ogr trying to tell me?  Am I using it properly, particularly
the +units flag?  Do I need some kind of explicit scale factor?

Q: Is there an easier way to manipulate / generate .prj's and .wkt's
aside from g.proj (so I don't necessarily have to start GRASS)?

Q: What's the preferred method of looking up EPSGs?  I rely on QGIS.

Q: Isn't this a generally useful and popular projection?  "Default
values for φ1 and φ2 are respectively 29.5 N and 45.5 N (values normally
used for maps of the conterminous United States)." (4)  Shouldn't it
have an easy-to-use EPSG?  Maybe it does and I just don't know where to
look.

Anyone still reading?  I will post again if I discover anything further,
but I have been obsessing about this for ~24 hours now, so I think I
will try something else.  Much obliged for any input from the group and
happy to provide omitted details.

Neil

(1)http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html
(2)http://www.intevation.de/rt/webrt?serial_num=3526&display=History
(3)http://www.remotesensing.org/geotiff/proj_list/guid7.html
(4) p. 35, Cartographic Projection Procedures for the UNIX Environment—A
    User’s Manual ftp://ftp.remotesensing.org/proj/new_docs/OF90-284.pdf

-------------- next part --------------
PROJCS["Albers Equal Area",
    GEOGCS["wgs84",
        DATUM["WGS_1984",
            SPHEROID["wgs84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["mi",1]]
-------------- next part --------------
INFO: Open of `temp'
using driver `ESRI Shapefile' successful.

Layer name: counties
Geometry: Polygon
Feature Count: 3140
Extent: (-6293474.061879, 2670863.569639) - (2256278.441196, 8557825.854853)
Layer SRS WKT:
PROJCS["Albers Equal Area",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["wgs84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["mi",1]]
NAME: String (32.0)
STATE_NAME: String (25.0)
STATE_FIPS: String (2.0)
CNTY_FIPS: String (3.0)
FIPS: String (5.0)
AREA: Real (18.4)
POP1990: Integer (10.0)
POP1997: Integer (10.0)
POP90_SQMI: Integer (6.0)
HOUSEHOLDS: Integer (8.0)
MALES: Integer (8.0)
FEMALES: Integer (8.0)
WHITE: Integer (8.0)
BLACK: Integer (8.0)
AMERI_ES: Integer (6.0)
ASIAN_PI: Integer (7.0)
OTHER: Integer (8.0)
HISPANIC: Integer (8.0)
AGE_UNDER5: Integer (7.0)
AGE_5_17: Integer (8.0)
AGE_18_29: Integer (8.0)
AGE_30_49: Integer (8.0)
AGE_50_64: Integer (8.0)
AGE_65_UP: Integer (7.0)
NEVERMARRY: Integer (8.0)
MARRIED: Integer (8.0)
SEPARATED: Integer (7.0)
WIDOWED: Integer (7.0)
DIVORCED: Integer (7.0)
HSEHLD_1_M: Integer (7.0)
HSEHLD_1_F: Integer (7.0)
MARHH_CHD: Integer (7.0)
MARHH_NO_C: Integer (7.0)
MHH_CHILD: Integer (6.0)
FHH_CHILD: Integer (7.0)
HSE_UNITS: Integer (8.0)
VACANT: Integer (7.0)
OWNER_OCC: Integer (8.0)
RENTER_OCC: Integer (8.0)
MEDIAN_VAL: Integer (7.0)
MEDIANRENT: Integer (4.0)
UNITS_1DET: Integer (8.0)
UNITS_1ATT: Integer (7.0)
UNITS2: Integer (7.0)
UNITS3_9: Integer (7.0)
UNITS10_49: Integer (7.0)
UNITS50_UP: Integer (7.0)
MOBILEHOME: Integer (6.0)
NO_FARMS87: Integer (5.0)
AVG_SIZE87: Integer (6.0)
CROP_ACR87: Integer (8.0)
AVG_SALE87: Integer (7.0)