[PROJ] About unit conversion

Kristian Evers kreve at sdfe.dk
Wed May 20 02:00:50 PDT 2020


I see that you don't specify the input units in your unitconvert steps. Try adding a
+xy_in parameter to your unitconvert steps and see if that fixes the problem.
When you only specify of of either +xy_in or +xy_out the operation becomes
somewhat uni-directional which is why not all of your pipelines works in reverse.

See https://proj.org/operations/conversions/unitconvert.html for more.

/Kristian

> -----Original Message-----
> From: PROJ <proj-bounces at lists.osgeo.org> On Behalf Of José Luis García
> Pallero
> Sent: 19. maj 2020 10:59
> To: proj at lists.osgeo.org
> Subject: Re: [PROJ] About unit conversion
> 
> 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(d
> ouble),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(doubl
> e),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!
> *****************************************


More information about the PROJ mailing list