proj=eqc with change of latitute

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

proj=eqc with change of latitute

Espen Isaksen
Hi!

This question originates from this earlier thread:

http://lists.umn.edu/cgi-bin/wa?A2=ind0801&L=mapserver-users&D=1&T=0&P=37923

But since I learned a lot during that thread and because the thread is
full of all kinds of questions now I would like to ask a very specific
question again in a new thread.

I have also checked out the archives at this list and the PROJ list. I
have found one similar question on this list in 2005:

http://lists.umn.edu/cgi-bin/wa?A2=ind0504&L=MAPSERVER-USERS&P=R75005

Here is my question:

I have set up a mapfile with projection parameters like this:

EXTENT 10.6661 59.9155 10.6824 59.9213
PROJECTION
    "init:epsg:4326"
END

From this mapfile I get a map over an small area in Oslo, Norway(I
cannot publish the data since the data is not free)

Then I change my map file to this:

EXTENT 37074.96 6669762.95 37982.21 6670408.60
PROJECTION
    "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +no_defs +units=m"
END

The extent is calculated by using PROJ like this:
proj "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +no_defs +units=m" < test.txt

where test.txt consists of the geographic coordinates in my first map file.

This mapfile creates a map about 18000 meters north of my first
mapfile. This is similar to the post from 2005. Longitude seems to be
correct.

I do not understand why Mapserver does not create approximatly the
same map in both map files? Is it not possible to use PROJ to convert
coordinates like this? If I set lat_ts=0 I do of course get the same
map, but that is not the projection I would like to use.

And if I cannot do it like this, how can I calculate the shift in the
y-coordinate(about 18000 meters in this example)?

Kind regards,
Espen
Reply | Threaded
Open this post in threaded view
|

Re: proj=eqc with change of latitute

Ed McNierney-4
Espen -

Those EXTENT values look odd.  Your projection's X-origin is still at 0 degrees longitude.  If your X coordinates are 10 degrees East, then the projected coordinates (in meters) shouldn't be about 37,000.  That would mean you're about 37 km east of the prime meridian, and you're not.  At 60 degrees North latitude one degree of longitude is 55.8 km, so if you're 10 degrees East you should have an X coordinate of about 558,000 meters (a little more, since you're more than 10 degrees East).

When I use PROJ or CS2CS to project your values I get exactly the same Y values but very different X values:

cs2cs +from +init=epsg:4326 +to +proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs +units=m
10.6661 59.9155
593672.41       6669762.95 0.00
10.6824 59.9213
594579.66       6670408.60 0.00

These X values seem reasonable for your projection.  The command line you give for PROJ shouldn't actually work as given - the quotes get in the way.  You should be using a command like the one above to CS2CS, or:

proj +proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs +units=m

for the conversion.  I'm a little puzzled as to how you're getting those results.  Could you try both the PROJ and CS2CS command lines above to double-check, and report the version of PROJ you're using?  There's probably something obvious here but I'm not seeing it.

     - Ed

Ed McNierney
Chief Mapmaker
Demand Media / TopoZone.com
73 Princeton Street, Suite 305
North Chelmsford, MA  01863
[hidden email]
Phone: +1 (978) 251-4242
Fax: +1 (978) 251-1396

-----Original Message-----
From: UMN MapServer Users List [mailto:[hidden email]] On Behalf Of Espen Isaksen
Sent: Tuesday, January 22, 2008 1:45 PM
To: [hidden email]
Subject: [UMN_MAPSERVER-USERS] proj=eqc with change of latitute

Hi!

This question originates from this earlier thread:

http://lists.umn.edu/cgi-bin/wa?A2=ind0801&L=mapserver-users&D=1&T=0&P=37923

But since I learned a lot during that thread and because the thread is
full of all kinds of questions now I would like to ask a very specific
question again in a new thread.

I have also checked out the archives at this list and the PROJ list. I
have found one similar question on this list in 2005:

http://lists.umn.edu/cgi-bin/wa?A2=ind0504&L=MAPSERVER-USERS&P=R75005

Here is my question:

I have set up a mapfile with projection parameters like this:

EXTENT 10.6661 59.9155 10.6824 59.9213
PROJECTION
    "init:epsg:4326"
END

From this mapfile I get a map over an small area in Oslo, Norway(I
cannot publish the data since the data is not free)

Then I change my map file to this:

EXTENT 37074.96 6669762.95 37982.21 6670408.60
PROJECTION
    "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +no_defs +units=m"
END

The extent is calculated by using PROJ like this:
proj "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
+datum=WGS84 +no_defs +units=m" < test.txt

where test.txt consists of the geographic coordinates in my first map file.

This mapfile creates a map about 18000 meters north of my first
mapfile. This is similar to the post from 2005. Longitude seems to be
correct.

I do not understand why Mapserver does not create approximatly the
same map in both map files? Is it not possible to use PROJ to convert
coordinates like this? If I set lat_ts=0 I do of course get the same
map, but that is not the projection I would like to use.

And if I cannot do it like this, how can I calculate the shift in the
y-coordinate(about 18000 meters in this example)?

Kind regards,
Espen
Reply | Threaded
Open this post in threaded view
|

Re: proj=eqc with change of latitute

Espen Isaksen
You are right, it was quite obvious. First of all, I made an error in
my e-mail since I was actually using lon_0=10 and not 0. But my main
problem was that proj was returning a wrong y-coordinate. The error is
the same as I mentioned in my first e-mail. Using cs2cs instead gave
me the correct coordinate, and those coordinates makes my map correct
:-)

Thanks for all you help Ed!

Espen

On 22/01/2008, Ed McNierney <[hidden email]> wrote:

> Espen -
>
> Those EXTENT values look odd.  Your projection's X-origin is still at 0 degrees longitude.  If your X coordinates are 10 degrees East, then the projected coordinates (in meters) shouldn't be about 37,000.  That would mean you're about 37 km east of the prime meridian, and you're not.  At 60 degrees North latitude one degree of longitude is 55.8 km, so if you're 10 degrees East you should have an X coordinate of about 558,000 meters (a little more, since you're more than 10 degrees East).
>
> When I use PROJ or CS2CS to project your values I get exactly the same Y values but very different X values:
>
> cs2cs +from +init=epsg:4326 +to +proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs +units=m
> 10.6661 59.9155
> 593672.41       6669762.95 0.00
> 10.6824 59.9213
> 594579.66       6670408.60 0.00
>
> These X values seem reasonable for your projection.  The command line you give for PROJ shouldn't actually work as given - the quotes get in the way.  You should be using a command like the one above to CS2CS, or:
>
> proj +proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +no_defs +units=m
>
> for the conversion.  I'm a little puzzled as to how you're getting those results.  Could you try both the PROJ and CS2CS command lines above to double-check, and report the version of PROJ you're using?  There's probably something obvious here but I'm not seeing it.
>
>      - Ed
>
> Ed McNierney
> Chief Mapmaker
> Demand Media / TopoZone.com
> 73 Princeton Street, Suite 305
> North Chelmsford, MA 01863
> [hidden email]
> Phone: +1 (978) 251-4242
> Fax: +1 (978) 251-1396
>
> -----Original Message-----
> From: UMN MapServer Users List [mailto:[hidden email]] On Behalf Of Espen Isaksen
> Sent: Tuesday, January 22, 2008 1:45 PM
> To: [hidden email]
> Subject: [UMN_MAPSERVER-USERS] proj=eqc with change of latitute
>
> Hi!
>
> This question originates from this earlier thread:
>
> http://lists.umn.edu/cgi-bin/wa?A2=ind0801&L=mapserver-users&D=1&T=0&P=37923
>
> But since I learned a lot during that thread and because the thread is
> full of all kinds of questions now I would like to ask a very specific
> question again in a new thread.
>
> I have also checked out the archives at this list and the PROJ list. I
> have found one similar question on this list in 2005:
>
> http://lists.umn.edu/cgi-bin/wa?A2=ind0504&L=MAPSERVER-USERS&P=R75005
>
> Here is my question:
>
> I have set up a mapfile with projection parameters like this:
>
> EXTENT 10.6661 59.9155 10.6824 59.9213
> PROJECTION
>     "init:epsg:4326"
> END
>
> From this mapfile I get a map over an small area in Oslo, Norway(I
> cannot publish the data since the data is not free)
>
> Then I change my map file to this:
>
> EXTENT 37074.96 6669762.95 37982.21 6670408.60
> PROJECTION
>     "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
> +datum=WGS84 +no_defs +units=m"
> END
>
> The extent is calculated by using PROJ like this:
> proj "+proj=eqc +lat_ts=60 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
> +datum=WGS84 +no_defs +units=m" < test.txt
>
> where test.txt consists of the geographic coordinates in my first map file.
>
> This mapfile creates a map about 18000 meters north of my first
> mapfile. This is similar to the post from 2005. Longitude seems to be
> correct.
>
> I do not understand why Mapserver does not create approximatly the
> same map in both map files? Is it not possible to use PROJ to convert
> coordinates like this? If I set lat_ts=0 I do of course get the same
> map, but that is not the projection I would like to use.
>
> And if I cannot do it like this, how can I calculate the shift in the
> y-coordinate(about 18000 meters in this example)?
>
> Kind regards,
> Espen
>