[PROJ] About unit conversion
José Luis García Pallero
jgpallero at gmail.com
Tue May 19 01:58:59 PDT 2020
I've written a small program in order to illustrate my last mail.
Attached (and pasted at the bottom of the message) I send the code.
The results are
Original data -> lon: -6.00000, lat: 43.00000 (degrees) -> lon:
-0.1047197551, lat: 0.7504915784 (radians)
rad to m -> X= 255466.981, Y=4765182.933 (+proj=utm +lon_0=3w +ellps=GRS80)
rad to km -> X= 255.467, Y= 4765.183 (+proj=pipeline +step
+proj=utm +lon_0=3w +ellps=GRS80 +step +proj=unitconvert +xy_out=km)
deg to m -> X= 255466.981, Y=4765182.933 (+proj=pipeline +step
+proj=unitconvert +xy_in=deg +step +proj=utm +lon_0=3w +ellps=GRS80)
deg to km -> X= 255.467, Y= 4765.183 (+proj=pipeline +step
+proj=unitconvert +xy_in=deg +step +proj=utm +lon_0=3w +ellps=GRS80
+step +proj=unitconvert +xy_out=km)
m to rad -> lon=-0.1047197551, lat= 0.7504915784 (+proj=utm +lon_0=3w
+ellps=GRS80)
km to rad -> lon=-0.0001306633, lat= 0.0000007501 (+proj=pipeline
+step +proj=unitconvert +xy_in=km +step +proj=utm +lon_0=3w
+ellps=GRS80)
m to deg -> lon=-0.1300126568, lat= 0.0130928532 (+proj=pipeline
+step +proj=utm +lon_0=3w +ellps=GRS80 +step +proj=unitconvert
+xy_out=deg)
km to deg -> lon=-0.0001307025, lat= 0.0000000131 (+proj=pipeline
+step +proj=unitconvert +xy_in=km +step +proj=utm +lon_0=3w
+ellps=GRS80 +step +proj=unitconvert +xy_out=deg)
As you can see, the only right INV conversion is the meters_to_radians
one. All other combinations produce all different results, but all
combinations for the FWD step are right. Im using PROJ 7.0.1 from the
Debian Sid repositories
***********************
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<proj.h>
#define CONST_PI (3.14159265358979323846264338328)
int main()
{
double lond1=-6.0,latd1=43.0,lond2=lond1,latd2=latd1;
double lonr1=lond1*CONST_PI/180.0,latr1=latd1*CONST_PI/180.0,lonr2=lonr1,latr2=latr1;
char param_rad_to_m[]="+proj=utm +lon_0=3w +ellps=GRS80";
char param_rad_to_km[]="+proj=pipeline +step +proj=utm +lon_0=3w
+ellps=GRS80 +step +proj=unitconvert +xy_out=km";
char param_deg_to_m[]="+proj=pipeline +step +proj=unitconvert
+xy_in=deg +step +proj=utm +lon_0=3w +ellps=GRS80";
char param_deg_to_km[]="+proj=pipeline +step +proj=unitconvert
+xy_in=deg +step +proj=utm +lon_0=3w +ellps=GRS80 +step
+proj=unitconvert +xy_out=km";
char param_m_to_rad[]="+proj=utm +lon_0=3w +ellps=GRS80";
char param_km_to_rad[]="+proj=pipeline +step +proj=unitconvert
+xy_in=km +step +proj=utm +lon_0=3w +ellps=GRS80";
char param_m_to_deg[]="+proj=pipeline +step +proj=utm +lon_0=3w
+ellps=GRS80 +step +proj=unitconvert +xy_out=deg";
char param_km_to_deg[]="+proj=pipeline +step +proj=unitconvert
+xy_in=km +step +proj=utm +lon_0=3w +ellps=GRS80 +step
+proj=unitconvert +xy_out=deg";
PJ_CONTEXT *C=proj_context_create();
PJ *proj=NULL;
size_t ld=sizeof(double);
//Display original data
printf("Original data -> lon: %.5lf, lat: %.5lf (degrees) -> lon:
%.10lf, lat: %.10lf (radians)\n\n",
lond1,latd1,lonr1,latr1);
//FWD, radiands to meters
proj = proj_create(C,param_rad_to_m);
proj_trans_generic(proj,PJ_FWD,&lonr1,ld,1,&latr1,ld,1,NULL,0,0,NULL,0,0);
printf("rad to m -> X=%11.3lf, Y=%11.3lf
(%s)\n",lonr1,latr1,param_rad_to_m);
proj_destroy(proj);
//FWD, radiands to kilometers
proj = proj_create(C,param_rad_to_km);
proj_trans_generic(proj,PJ_FWD,&lonr2,ld,1,&latr2,ld,1,NULL,0,0,NULL,0,0);
printf("rad to km -> X=%11.3lf, Y=%11.3lf
(%s)\n",lonr2,latr2,param_rad_to_km);
proj_destroy(proj);
//FWD, degrees to meters
proj = proj_create(C,param_deg_to_m);
proj_trans_generic(proj,PJ_FWD,&lond1,ld,1,&latd1,ld,1,NULL,0,0,NULL,0,0);
printf("deg to m -> X=%11.3lf, Y=%11.3lf
(%s)\n",lond1,latd1,param_deg_to_m);
proj_destroy(proj);
//FWD, degrees to kilometers
proj = proj_create(C,param_deg_to_km);
proj_trans_generic(proj,PJ_FWD,&lond2,ld,1,&latd2,ld,1,NULL,0,0,NULL,0,0);
printf("deg to km -> X=%11.3lf, Y=%11.3lf
(%s)\n\n",lond2,latd2,param_deg_to_km);
proj_destroy(proj);
//INV, meters to radians
proj = proj_create(C,param_m_to_rad);
proj_trans_generic(proj,PJ_INV,&lonr1,ld,1,&latr1,ld,1,NULL,0,0,NULL,0,0);
printf("m to rad -> lon=%13.10lf, lat=%13.10lf
(%s)\n",lonr1,latr1,param_m_to_rad);
proj_destroy(proj);
//INV, kilometers to radians
proj = proj_create(C,param_km_to_rad);
proj_trans_generic(proj,PJ_INV,&lonr2,ld,1,&latr2,ld,1,NULL,0,0,NULL,0,0);
printf("km to rad -> lon=%13.10lf, lat=%13.10lf
(%s)\n",lonr2,latr2,param_km_to_rad);
proj_destroy(proj);
//INV, meters to degrees
proj = proj_create(C,param_m_to_deg);
proj_trans_generic(proj,PJ_INV,&lond1,ld,1,&latd1,ld,1,NULL,0,0,NULL,0,0);
printf("m to deg -> lon=%13.10lf, lat=%13.10lf
(%s)\n",lond1,latd1,param_m_to_deg);
proj_destroy(proj);
//INV, kilometers to degrees
proj = proj_create(C,param_km_to_deg);
proj_trans_generic(proj,PJ_INV,&lond2,ld,1,&latd2,ld,1,NULL,0,0,NULL,0,0);
printf("km to deg -> lon=%13.10lf, lat=%13.10lf
(%s)\n",lond2,latd2,param_km_to_deg);
proj_destroy(proj);
proj_context_destroy(C);
return 0;
}
El lun., 18 may. 2020 a las 23:03, José Luis García Pallero
(<jgpallero at gmail.com>) escribió:
>
> Hello:
>
> I'm trying to do dome projections (forward and inverse) using the
> function proj_trans_generic() and parameter string in PROJ format (+
> style) and I'm a bit confused about unit conversions. My reference for
> the syntax is https://proj.org/operations/conversions/unitconvert.html
>
> For the forward problem I write
>
> C = proj_context_create();
> projection = proj_create(C,param);
> proj_trans_generic(projection,PJ_FWD,lon,sizeof(double),nElem,lat,sizeof(double),nElem,NULL,0,0,NULL,0,0);
>
> I've used the param string
>
> '+proj=pipeline +step +proj=unitconvert +xy_in=deg +step +proj=utm
> +lon_0=3w +ellps=GRS80 +step +proj=unitconvert +xy_out=km'
>
> which means that the input coordinates are in degrees and the output
> in kilometers. All works right. The results are the same if I pass the
> input coordinates in radians and I use the param string
>
> '+proj=pipeline +step +proj=utm +lon_0=3w +ellps=GRS80 +step
> +proj=unitconvert +xy_out=km'
>
> The problem comes when I try to do the inverse step. The procedure is the same:
>
> C = proj_context_create();
> projection = proj_create(C,param);
> proj_trans_generic(projection,PJ_INV,x,sizeof(double),nElem,y,sizeof(double),nElem,NULL,0,0,NULL,0,0);
>
> In the input xoordinates are in meters and I use the param string
>
> '+proj=utm +lon_0=3w +ellps=GRS80'
>
> The resulting longtude and latitude are correct and the units are
> radians. But if I use the param string
>
> '+proj=pipeline +step +proj=utm +lon_0=3w +ellps=GRS80 +step
> +proj=unitconvert +xy_out=deg'
>
> The output is totally wrong. I expect degrees, but the values have no
> sense. Also, if I use the input X and Y in kilometers ans the param
> string is
>
> '+proj=pipeline +step +proj=unitconvert +xy_in=km +step +proj=utm
> +lon_0=3w +ellps=GRS80'
>
> the result is wrong again.
>
> I suppose I'm doing something wrong with unit conversion, but I can not find it
>
> Can someone help me?
>
> Thanks
>
> --
> *****************************************
> José Luis García Pallero
> jgpallero at gmail.com
> (o<
> / / \
> V_/_
> Use Debian GNU/Linux and enjoy!
> *****************************************
--
*****************************************
José Luis García Pallero
jgpallero at gmail.com
(o<
/ / \
V_/_
Use Debian GNU/Linux and enjoy!
*****************************************
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test-proj.c
Type: text/x-csrc
Size: 3604 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200519/918e3c04/attachment-0001.c>
More information about the PROJ
mailing list