[Proj] Shapefile and EXTENTS

Wed May 11 12:16:44 EDT 2005

Martin,

Thanks.  I will take a look at OGR.  This helps a whole lot!

Matt

On May 10, 2005, at 9:26 PM, Chapman, Martin wrote:

> 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
>
>
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/proj