Update:
When I work in GRASS and use g.setproj I can tell it that my units
should be miles in the prompts. The result is:
> g.proj -jf
+proj=aea +lat_0=23.0000000000 +lat_1=29.5000000000 +lat_2=45.5000000000
+lon_0=-96.0000000000 +x_0=0.0000000000 +y_0=0.0000000000 +a=6378137
+rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000
> g.proj -w
PROJCS["Albers Equal Area",
GEOGCS["wgs84",
DATUM["WGS_1984",
SPHEROID["wgs84",6378137,298.257223563],
TOWGS84[0.000,0.000,0.000]],
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",23],
PARAMETER["longitude_of_center",-96],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["miles",1609.344]]
> g.proj -fw > aeaUS.wkt
>ogr2ogr -s_srs EPSG:4326 -t_srs aeaUS.wkt -where 'state_fips not in
("02","15")' temp usagen counties
# leaving out AK and HI -- sorry folks!
QGIS doesn't like the proj4 string for a custom projection (it complains
about the lack of an +ellps[sic] clause), but it picks it up from
the .prj generated by ogr2ogr just fine. Now I should be able to
generate my buffers in GRASS.
Is there a bug somewhere among all of this? The key appears to be the
UNIT clause of the WKT. Comments?
g.setproj told me to do a g.region -d but it didn't seem to do the right
thing and I had to start over with a new location. I know this isn't a
GRASS list, but while I have the floor . . .
On Fri, 2006-06-16 at 10:22 -0500, Neil Best wrote:
> 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
>
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/gdal-dev