[Proj] Shapefile and EXTENTS

Tue May 10 22:26:24 EDT 2005

Matt,

Glad to help.  I was guessing that was your problem, I had the same issue once.  

To answer your question, projected spatial systems like LAE (Lambert Azimuthal Equal) are only in linear units.  The reason why, is because they were originally angular units (Geographic Coordinate System), and the process of "Projecting" the points translates them into a Projected Coordinate System (LAE in your case).  So, you can't have a projected system in angular units, but you can translate your points back and forth between linear and angular units, or in other words, between a "Geographic Coordinate System" and a "Projected Coordinate System".  This is called a coordinate transform operation.  I would suggest looking at http://www.remotesensing.org/ogr .  It's an open source project for vector data that has a higher level abstraction library for dealing with spatial reference systems that uses proj4 under the hood.  There are many useful classes and utilities for spatial operations, including a simple to use OGRCoordinateTransform object that will help transform point data from one coordinate space to another.  There are helpful code tutorials and the library is easy to download and install.  Also, you can read more about how coordinate systems work.  I hope this helps.  Let me know if you need any additional code examples.

One last thing.  To get your LAE points into angular units in Proj you would convert them to a Geographic system like WGS84 something like this where the sInitString variable is a proj string that identifies the LAE projection:

m_sInitString = sInitString;

char seps[] = "+ ";

char *token = NULL;

vector< string> vTokens;

char* pszInit = (char*) sInitString.c_str();

token = strtok( pszInit, seps );

while( token != NULL )

{

char buffer[255];

strcpy(buffer, token);

string str(buffer);

vTokens.push_back(str);

token = strtok( NULL, seps );

}

char* parms[40] = { new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255],

new char[255], new char[255], new char[255], new char[255], new char[255] };

for (int j = 0; j < 40; j++)

strncpy(parms[j], "\n", 255);

for (int i = 0; i < vTokens.size(); i++)

{

char* psz = (char*) vTokens[i].c_str();

strcpy(parms[i], psz);

}

if (m_pProjection != NULL)

pj_free(m_pProjection);

if ( !(m_pProjection = pj_init(sizeof(parms) / sizeof(char*), parms)) ) 

throw CString("Projection initialization failed");

for (int j = 0; j < 40; j++)

delete parms[j];

// then for each point call inverse

vector< float> vPoint;

projUV point;

point.v = (double) nY;

point.u = (double) nX;

point = pj_inv(point, m_pProjection);

I think this is how you do it in proj4.  If that doesn't work then call pj_fwd instead, but I'm pretty sure that's it.  It's hard to remember since I now use OGR exclusively for my spatial reference transforms.

Martin

________________________________

From: proj-bounces at xserve.flids.com on behalf of Matthew Roberson
Sent: Tue 5/10/2005 8:00 AM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Shapefile and EXTENTS

Martin,

Thank you.  After reading your reply, I understand why I can't see my
map. Another question I have is, can I create a Lambert Azimuthal Equal
Area projection in Angular units using the proj utility? How would I do
that?  Thanks

Matt

On May 9, 2005, at 5:16 PM, Chapman, Martin wrote:

> Matt,
>
> The issue is that you are converting from angular units to linear units
> when projecting from geographic to a projected system.  Angular units
> range from -180 to 180 and 90 to -90 degrees.  After the re-projection
> your values are in false easting and false northing which are huge
> number offsets to put all values in the positive quadrant of the
> cartesian plane.  If you are trying to display linear units in a map
> configured to show angular units, then your linear units will be way
> off
> the screen.  Make sure your map is rendering in the same spatial
> reference system as the file you are trying to display, or vice-versa.
>
> Martin
>
> -----Original Message-----
> From: proj-bounces at xserve.flids.com
> [mailto:proj-bounces at xserve.flids.com] On Behalf Of Matthew Roberson
> Sent: Monday, May 09, 2005 3:43 PM
> To: proj at xserve.flids.com
> Subject: [Proj] Shapefile and EXTENTS
>
>
> I am trying to do a project a map of the state of Missouri that was
> created by the US Census Bureau(zt29_d00.shp).
>
> I used ogrinfo to get the following information about the EXTENT:
>
> [imac-dev03:~/Documents] matthewr% ogrinfo -so zt29_d00.shp zt29_d00
> INFO: Open of `zt29_d00.shp'
> using driver `ESRI Shapefile' successful.
>
> Layer name: zt29_d00
> Geometry: Polygon
> Feature Count: 1341
> Extent: (-95.774704, 35.995683) - (-89.098843, 40.613640)
> Layer SRS WKT:
> (unknown)
> AREA: Real (20.5)
> PERIMETER: Real (20.5)
> ZT29_D00_: Integer (11.0)
> ZT29_D00_I: Integer (11.0)
> ZCTA: String (5.0)
> NAME: String (90.0)
> LSAD: String (2.0)
> LSAD_TRANS: String (50.0)
>
> Then, I used proj4 to find lat long values for the a Lambert Azimuthal
> Equal Area projection.
>
> [imac-dev03:~/Documents] matthewr% proj +proj=laea +ellps=clrk66 +R_A
> -95.774704 35.995683
> -7567149.62     5525023.18
> -89.098843 40.613640
> -6798293.00     5830370.88
>
> When I use these values with MapServer to make a map, nothing shows up.
>   Am I missing something.  Should I post this question in this list or
> in the MapServer User's group? Any help would be appreciated.
>
> Thanks,
>
> Matt
>
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com http://xserve.flids.com/mailman/listinfo/proj
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/proj

_______________________________________________
Proj mailing list
Proj at xserve.flids.com
http://xserve.flids.com/mailman/listinfo/proj