Projection conversion problem

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

Projection conversion problem

Damian Rakowski
Hi,

my names is Damian, I'm new to proj and projections topic.
I need to perform transformation between two projections, and I have a problem defining projection string properly.

I have coordinates in WGS84 like:

2.2945 48.858222

and I need to transform them to WGS84-like projection with different unit - encoded on integer.
Basically what I need to do is to transform coordinates in the following way:

floor(coordinate * 2^32/360)
2^32/360 = 11930464,711111111111111111111111


After conversion I expect to get the following coordinates:

27374451 582901293


I tried:

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs  +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +k_0=11930464,711111111111111111111111 << EOF
> 2.2945 48.858222
> EOF
2d17'40.2"E    48d51'29.599"N 0.000


cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs  +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +k=11930464,711111111111111111111111 << EOF
>2.2945 48.858222
>EOF

2d17'40.2"E    48d51'29.599"N 0.000

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs  +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +m=11930464,711111111111111111111111 << EOF
2.2945 48.858222
EOF

The same for reverse conversion:

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +k_0=11930464,711111111111111111111111  +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs << EOF        
27374451 582901293
EOF

27374451dE    582901293dN 0.000


So it seems that k, k_0, m does not bring any effect. I must not understand or missing something. Could you please help me?

Best Regards,

Damian

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

peifer
On 2016-06-27 11:43, Damian Rakowski wrote:

>
> I have coordinates in WGS84 like:
>
> 2.2945 48.858222
>
> and I need to transform them to WGS84-like projection with different
> unit - encoded on integer.
> Basically what I need to do is to transform coordinates in the following
> way:
>
> floor(coordinate * 2^32/360)
> 2^32/360 = 11930464,711111111111111111111111
>
>
> After conversion I expect to get the following coordinates:
>
> 27374451 582901293
>

This looks like Garmin's "semicircle units". You don't need proj for
calculating them. Any of awk/perl/python/ruby/... will do.

Hermann
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Andre Joost
In reply to this post by Damian Rakowski
Am 27.06.2016 um 11:43 schrieb Damian Rakowski:

>
> So it seems that k, k_0, m does not bring any effect. I must not
> understand or missing something. Could you please help me?
>

These scale factors are only useful for certain projected coordinate
systems, like tmerc.

Since you are projecting from longlat to longlat, they will not be used
in the underlying code.

Just scale your coordinates without using proj.

HTH,
André Joost


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Damian Rakowski
In reply to this post by peifer
Thank you for your reply. That's true i can convert it using scripts,
but I was thinking about different use case. I had a hope that thanks to
proj integration in postgis I could define new coordinate system with
proper proj configuration string (similar to this one:
http://spatialreference.org/ref/epsg/4326/postgis/), and use such
encoded coordinates directly in postgis.

Best Regards,

Damian

W dniu 2016-06-27 o 20:14, Hermann Peifer pisze:

> On 2016-06-27 11:43, Damian Rakowski wrote:
>> I have coordinates in WGS84 like:
>>
>> 2.2945 48.858222
>>
>> and I need to transform them to WGS84-like projection with different
>> unit - encoded on integer.
>> Basically what I need to do is to transform coordinates in the following
>> way:
>>
>> floor(coordinate * 2^32/360)
>> 2^32/360 = 11930464,711111111111111111111111
>>
>>
>> After conversion I expect to get the following coordinates:
>>
>> 27374451 582901293
>>
> This looks like Garmin's "semicircle units". You don't need proj for
> calculating them. Any of awk/perl/python/ruby/... will do.
>
> Hermann

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Norman Vine
I believe you are looking for the plate caree projection
https://epsg.io/32662

   from: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
   to:     +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

"eqc" above is "Equidistant Cylindrical (Plate Caree)" - which is an alias of OGC WKT Name: "Equirectangular"

You will need to truncate this to integer

Norman


> On Jun 28, 2016, at 11:43 AM, Damian Rakowski <[hidden email]> wrote:
>
> Thank you for your reply. That's true i can convert it using scripts,
> but I was thinking about different use case. I had a hope that thanks to
> proj integration in postgis I could define new coordinate system with
> proper proj configuration string (similar to this one:
> http://spatialreference.org/ref/epsg/4326/postgis/), and use such
> encoded coordinates directly in postgis.
>
> Best Regards,
>
> Damian
>
> W dniu 2016-06-27 o 20:14, Hermann Peifer pisze:
>> On 2016-06-27 11:43, Damian Rakowski wrote:
>>> I have coordinates in WGS84 like:
>>>
>>> 2.2945 48.858222
>>>
>>> and I need to transform them to WGS84-like projection with different
>>> unit - encoded on integer.
>>> Basically what I need to do is to transform coordinates in the following
>>> way:
>>>
>>> floor(coordinate * 2^32/360)
>>> 2^32/360 = 11930464,711111111111111111111111
>>>
>>>
>>> After conversion I expect to get the following coordinates:
>>>
>>> 27374451 582901293
>>>
>> This looks like Garmin's "semicircle units". You don't need proj for
>> calculating them. Any of awk/perl/python/ruby/... will do.
>>
>> Hermann
>
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

peifer
On 2016-06-28 18:06, Norman Vine wrote:

> I believe you are looking for the plate caree projection
> https://epsg.io/32662
>
>    from: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
>    to:     +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
>
> "eqc" above is "Equidistant Cylindrical (Plate Caree)" - which is an alias of OGC WKT Name: "Equirectangular"
>
> You will need to truncate this to integer
>

With the given sample coordinates, this will give:

$ echo 2.2945 48.858222 | cs2cs -v +proj=longlat +datum=WGS84 +no_defs
+to +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs
# ---- From Coordinate System ----
#Lat/long (Geodetic alias)
#
# +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# ---- To Coordinate System ----
#Equidistant Cylindrical (Plate Caree)
# Cyl, Sph
# lat_ts=[, lat_0=0]
# +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m
# +no_defs +ellps=WGS84 +towgs84=0,0,0
255422.57 5438872.39 0.00

The expected result was: 27374451 582901293

Hermann
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Andre Joost
Am 28.06.2016 um 19:40 schrieb Hermann Peifer:

>
> With the given sample coordinates, this will give:
>
> $ echo 2.2945 48.858222 | cs2cs -v +proj=longlat +datum=WGS84 +no_defs
> +to +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs
> # ---- From Coordinate System ----
> #Lat/long (Geodetic alias)
> #
> # +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
> # ---- To Coordinate System ----
> #Equidistant Cylindrical (Plate Caree)
> # Cyl, Sph
> # lat_ts=[, lat_0=0]
> # +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> # +no_defs +ellps=WGS84 +towgs84=0,0,0
> 255422.57 5438872.39 0.00
>
> The expected result was: 27374451 582901293
>

I suggest to take a sphere, and scale with the units:

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000
+to_meter=0.0093202511 +no_defs

which results to 27374451.23 582901292.36 0.00

HTH,
André Joost


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Damian Rakowski
That  looks very interesting, thank you for your reply.

One more question:

6371000 - this is an Earth radius
0.0093202511 - how did you got this value? I just want to understand :)

Best Regards,

Damian


19:18 środa, 2016-6-29, Andre Joost <[hidden email]> napisał(a):


Am 28.06.2016 um 19:40 schrieb Hermann Peifer:

>
> With the given sample coordinates, this will give:
>
> $ echo 2.2945 48.858222 | cs2cs -v +proj=longlat +datum=WGS84 +no_defs
> +to +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs
> # ---- From Coordinate System ----
> #Lat/long (Geodetic alias)
> #
> # +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
> # ---- To Coordinate System ----
> #Equidistant Cylindrical (Plate Caree)
> #    Cyl, Sph
> #    lat_ts=[, lat_0=0]
> # +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> # +no_defs +ellps=WGS84 +towgs84=0,0,0
> 255422.57    5438872.39 0.00
>
> The expected result was: 27374451 582901293
>

I suggest to take a sphere, and scale with the units:

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000
+to_meter=0.0093202511 +no_defs

which results to 27374451.23    582901292.36 0.00

HTH,
André Joost



_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj



_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Andre Joost
Am 30.06.2016 um 10:20 schrieb Damian Rakowski:
> That  looks very interesting, thank you for your reply. One more
> question: 6371000 - this is an Earth radius0.0093202511 - how did you
> got this value? I just want to understand :) Best Regards, Damian
>


In a first round, I substituted the ellipsoid by the radius:

echo 2.2945 48.858222   | cs2cs -v +proj=longlat +R=6371000 +no_defs +to
+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m
+no_defs

255136.76 5432786.41

then divided the result in meters by the expected result. Since it is
the same for both coordinates, it seems to be right.

 From John P. Snyders "Map Projections - A Working Manual" p.91, you can
take the exact formulas for the sphere:

x = R (lambda-lambda0) cos phi1
y = R phi

with angles in radians. This is linear to the angles, and you just have
to compare them with your scale factor.


Feel free to test it with other coordinate pairs around the globe.

HTH,
André Joost

> Am 28.06.2016 um 19:40 schrieb Hermann Peifer:
>
>>
>> With the given sample coordinates, this will give:
>>
>> $ echo 2.2945 48.858222 | cs2cs -v +proj=longlat +datum=WGS84
>> +no_defs +to +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0
>> +datum=WGS84 +units=m +no_defs # ---- From Coordinate System ----
>> #Lat/long (Geodetic alias) # # +proj=longlat +datum=WGS84 +no_defs
>> +ellps=WGS84 +towgs84=0,0,0 # ---- To Coordinate System ----
>> #Equidistant Cylindrical (Plate Caree) #    Cyl, Sph #    lat_ts=[,
>> lat_0=0] # +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0
>> +datum=WGS84 +units=m # +no_defs +ellps=WGS84 +towgs84=0,0,0
>> 255422.57    5438872.39 0.00
>>
>> The expected result was: 27374451 582901293
>>
>
> I suggest to take a sphere, and scale with the units:
>
> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000
> +to_meter=0.0093202511 +no_defs
>
> which results to 27374451.23    582901292.36 0.00
>
> HTH, André Joost
>
>
> _______________________________________________ Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>
>
>
>
>
> _______________________________________________ Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection conversion problem

Andre Joost
Am 30.06.2016 um 13:24 schrieb Andre Joost:
> Am 30.06.2016 um 10:20 schrieb Damian Rakowski:
>> That  looks very interesting, thank you for your reply. One more
>> question: 6371000 - this is an Earth radius0.0093202511 - how did you
>> got this value? I just want to understand :) Best Regards, Damian
>>

>
>   From John P. Snyders "Map Projections - A Working Manual" p.91, you can
> take the exact formulas for the sphere:
>
> x = R (lambda-lambda0) cos phi1
> y = R phi
>
> with angles in radians. This is linear to the angles, and you just have
> to compare them with your scale factor.
>

So my to_meter factor is:

PI() * 6371000 / 2^31 = 0.009320251

Greetings,
André Joost


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj