---

gtkfmapsproj.c

Go to the documentation of this file.

/* gtkfmapsproj.c - Projection library of gtkfmaps
* 
 *  Copyright (c) 2000, Franck Martin <franck@sopac.org>, Eric G. Miller <egm2@jps.net>
 *
 * This Program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This Program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this Program; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*
 * Modified by the FMaps Team and others 2000.  See the AUTHORS
 * file for a list of people on the FMaps Team.  See the ChangeLog
 * files for a list of changes.  These files are distributed with
 * FMaps at http://FMaps.sourceforge.net/.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <float.h>
#include <string.h>

#include "gtkfmaps.h"
#include "gtkfmapsproj.h"


/*
 *  Delimiters for input and output strings.  LDELIM, RDELIM, and DELIM are
 *  left, right, and separator delimiters, respectively.  LDELIM_EP,
 *  RDELIM_EP are left and right delimiters for paths with endpoints.
 *
 */

#define LDELIM                  '('
#define RDELIM                  ')'
#define DELIM                   ','
#define LDELIM_EP               '['
#define RDELIM_EP               ']'
#define LDELIM_C                '<'
#define RDELIM_C                '>'


static int point_decode(char *str, gdouble *x, gdouble *y, gdouble *z, char **s);

void
gtk_fmaps_display_engine (GtkFMaps *fmaps)
{
  gint i;
  gint j;
  gchar *buffer;
  PGresult *qryResult;
  gint nbtuples;
  Point *points;
  gint nbpoints;
  gfloat r,g,b,a;
  gfloat size;
  gint type=-1;

  g_return_if_fail (fmaps !=NULL);
  g_return_if_fail (GTK_IS_FMAPS (fmaps));

  for(i=0;i<fmaps->nbtables;i++)
  {
        if (fmaps->table[i]->visible)
        {
          if (fmaps->table[i]->vector)
          {
                    buffer=g_strdup_printf("SELECT %s from f_geo_%s ",F_GEO_TABLE_FIELD, fmaps->table[i]->name);
              qryResult = PQexec(fmaps->conn,buffer);

              /* debug */
              if (!(g_str_equal(PQresultErrorMessage(qryResult),"")))
              {
                 g_warning("%s : %s",buffer, PQresultErrorMessage(qryResult));
              }
        
                  g_free(buffer);
                  nbtuples = PQntuples(qryResult);
                  j=0;
                          while(j<nbtuples)
                  {
                      buffer = PQgetvalue(qryResult,j,10);
                      points=gtk_fmaps_get_points(fmaps,i,buffer,&nbpoints,&type);      
                      size=strtod(PQgetvalue(qryResult,j,1),NULL);
                
                      r=strtod(PQgetvalue(qryResult,j,4),NULL);
                      g=strtod(PQgetvalue(qryResult,j,5),NULL);
                      b=strtod(PQgetvalue(qryResult,j,6),NULL);
                      a=strtod(PQgetvalue(qryResult,j,7),NULL);
        
                      if (nbpoints>0 && type!=-1)
                      {
                        glColor4f(r,g,b,a);
                        glBegin(type);
                        glPointSize(size);
                        for (i=0;i<nbpoints;i++)
                        {               
                          glVertex3f(points[i].x,points[i].y,points[i].z);
                        }
                        glEnd();
                      }
                      g_free(points);
                      j++;
                    }
                    PQclear(qryResult);         
                }
          }
        }
}

Point* gtk_fmaps_get_points(GtkFMaps *fmaps, gint tablenb, gchar *buffer, gint *nbpoints, gint *type)
{
  gint i;
  gdouble x;
  gdouble y;
  gdouble z;
  gint pointnb;
  gint npoints;
  Point *points=NULL;
  gboolean result=FALSE;
  gchar *str=NULL;

  g_return_val_if_fail(fmaps!=NULL,0);
  g_return_val_if_fail(buffer!=NULL,0);

  x=0;
  y=0;
  z=0;
  i=0;
  pointnb=0;
  npoints=0;

    
  if (strncmp("POINT",buffer,5)==0)
  {
    *type=GL_POINTS;
    i=5;
  } else if (strncmp("POLYLINE",buffer,8)==0)
  {
    *type=GL_LINE_STRIP;
    i=8;
  } else
  {
    *type=-1;
  }

  while(i<strlen(buffer))
  {
    if (buffer[i]==LDELIM)
    {
      str=g_strdup(buffer+i);
      result=point_decode(str,&x,&y,&z,NULL);
      g_free(str);
      npoints++;
      if (npoints==1)
      {
        points=g_new(Point,npoints);
      } else
      {
        points=g_renew(Point,points,npoints);
      }
      points[npoints-1].x=x;
      points[npoints-1].y=y;
      points[npoints-1].z=z;

   }

  i++;
  };

  *nbpoints=npoints;
  return(points);

}


/* Private Functions */
static int
point_decode(char *str, gdouble *x, gdouble *y, gdouble *z, char **s)
{
        int     has_delim;
        char    *cp;

        if(NULL == str)
                return FALSE;

        while(isspace(*str))
                str++;
        if ((has_delim = (*str == LDELIM)))
                str++;
        while(isspace(*str))
                str++;
        *x = strtod(str, &cp);
        if (cp <= str)
                return FALSE;
        while(isspace(*cp))
                cp++;
        if (*cp++ != DELIM)
                return FALSE;
        while (isspace(*cp))
                cp++;
        *y = strtod(cp, &str);
        if(str <= cp)
                return FALSE;
        while(isspace(*str))
                str++;
        if (*str++ != DELIM)
                return FALSE;
        while(isspace(*str))
                str++;
        *z = strtod(str, &cp);
        if(cp <= str)
                return FALSE;
        while(isspace(*cp))
                cp++;
        if(has_delim)
        {
                if (*cp != RDELIM)
                        return FALSE;
                cp++;
                while (isspace(*cp))
                        cp++;
        }
        if (s != NULL)
                *s = cp;

        return TRUE;
} /* point_decode() */

gint16 gtk_fmaps_displayx(gdouble x, GtkFMaps *fmaps)   
{
  gint16 xx;

  g_return_val_if_fail(fmaps!=NULL,0);

  xx=((gint16)((((x-fmaps->centerx)/fmaps->width)*GTK_WIDGET(fmaps)->allocation.width)+GTK_WIDGET(fmaps)->allocation.width/2));
  return(xx);
}

gint16 gtk_fmaps_displayy(gdouble y, GtkFMaps *fmaps)
{
  gint16 yy;
  g_return_val_if_fail(fmaps!=NULL,0);

  yy=((gint16)((((y-fmaps->centery)/fmaps->width)*GTK_WIDGET(fmaps)->allocation.width)+GTK_WIDGET(fmaps)->allocation.height/2));
  return(yy);
}

/* convert longitude, latitude and ellipsoid height into cartesian coordinates */
void gtk_fmaps_lonlath_2_xyz(GtkFMapsProjDatum *projdatum, gdouble lon, gdouble lat, gdouble h, gdouble *x, gdouble *y, gdouble *z)
{
  gdouble N;

  g_return_if_fail(projdatum != NULL);

  N=projdatum->a/sqrt(1-projdatum->e*projdatum->e*sin(lat)*sin(lat));
  *x=(N+h)*cos(lat)*cos(lon);
  *y=(N=h)*cos(lat)*sin(lon);
  *z=((1-projdatum->e*projdatum->e)*N+h)*sin(lat);

}

/* convert cartesian coordinates into longitude, latitude and ellipsoid height */
void gtk_fmaps_xyz_2_lonlath(GtkFMapsProjDatum *projdatum, gdouble x, gdouble y, gdouble z, gdouble *lon, gdouble *lat, gdouble *h)
{

  gdouble R;
  gdouble N;
  gdouble lat1;

  g_return_if_fail(projdatum!=NULL);

  R=sqrt(x*x+y*y);
  N=0;
  if (x==0)
  {
    *lon=atan(y/x);
  } else
  {
    if (y>=0)
    {
      *lon=M_PI_2;
    } else
    {
      *lon=-M_PI_2;
    }
  }

  if (x<0)
  {
    if (y<0)
    {
      *lon=*lon-M_PI;
    } else
    {
      *lon=*lon+M_PI;
    }
  }
  lat1=0;
  if (R!=0)
  {
    do
    {
      lat1=*lat;
      N=projdatum->a/sqrt(1-projdatum->e*projdatum->e*sin(lat1)*sin(lat1));
      *lat=atan((z+projdatum->e*projdatum->e*N*sin(lat1))/R);
      *h=R/cos(lat1)-N;
    } while (lat1-*lat!=0);
  } else
  {
    if (z>=0)
    {
      *lat=M_PI_2;
    } else
    {
      *lat=-M_PI_2;
    }
    *h=-N;
  }
}

/* transform cartesian coordinates from one datum to WGS84 datun if forward otherwise from WGS84 */
void gtk_fmaps_bursa_wolfe(GtkFMapsProjDatum *projdatum, gboolean forward, gdouble *x, gdouble *y, gdouble *z)
{
  gdouble x1,y1,z1;

  g_return_if_fail(projdatum!=NULL);

  x1=*x;
  y1=*y;
  z1=*z;

  if (forward)
  {
    *x=projdatum->r11*x1+projdatum->r12*y1+projdatum->r13*z1+projdatum->dx;
    *y=projdatum->r21*x1+projdatum->r22*y1+projdatum->r23*z1+projdatum->dy;
    *z=projdatum->r31*x1+projdatum->r32*y1+projdatum->r33*z1+projdatum->dz;
  } else
  {
    *x=projdatum->r11*x1+projdatum->r12*y1+projdatum->r13*z1-projdatum->dx;
    *y=projdatum->r21*x1+projdatum->r22*y1+projdatum->r23*z1-projdatum->dy;
    *z=projdatum->r31*x1+projdatum->r32*y1+projdatum->r33*z1-projdatum->dz;
  }
}

/* transform cartesian coordinates from one datum to WGS84 datun if forward otherwise from WGS84 */
void gtk_fmaps_molendensky(GtkFMapsProjDatum *projdatum, gboolean forward, gdouble *x, gdouble *y, gdouble *z)
{
  g_return_if_fail(projdatum!=NULL);

  if (forward)
  {
    *x=*x+projdatum->dx;
    *y=*y+projdatum->dy;
    *z=*z+projdatum->dz;
  } else
  {
    *x=*x-projdatum->dx;
    *y=*y-projdatum->dy;
    *z=*z-projdatum->dz;
  }
}

/* initialise the datum transformation parameters */
void gtk_fmaps_projection_init(GtkFMapsProjDatum *projdatum)
{
  g_return_if_fail(projdatum!=NULL);

  switch(projdatum->datumtransform)
  {
  case DAT_NONE:
  case DAT_MOLENDENSKY:
    projdatum->r11=1;
    projdatum->r12=0;
    projdatum->r13=0;

    projdatum->r21=0;
    projdatum->r22=1;
    projdatum->r23=0;

    projdatum->r31=0;
    projdatum->r32=0;
    projdatum->r33=1;
    break;
  case DAT_BURSA_WOLFE:
    projdatum->r11=projdatum->m*cos(projdatum->ez)*cos(projdatum->ey);
    projdatum->r21=-projdatum->m*sin(projdatum->ez)*cos(projdatum->ey);
    projdatum->r31=projdatum->m*sin(projdatum->ey);

    projdatum->r12=projdatum->m*(cos(projdatum->ez)*sin(projdatum->ey)*sin(projdatum->ex)+sin(projdatum->ez)*cos(projdatum->ex));
    projdatum->r22=projdatum->m*(cos(projdatum->ez)*cos(projdatum->ex)-sin(projdatum->ez)*sin(projdatum->ey)*sin(projdatum->ex));
    projdatum->r32=-projdatum->m*cos(projdatum->ey)*sin(projdatum->ex);

    projdatum->r13=projdatum->m*(sin(projdatum->ey)*sin(projdatum->ex)-cos(projdatum->ey)*sin(projdatum->ey)*cos(projdatum->ex));
    projdatum->r23=projdatum->m*(sin(projdatum->ez)*sin(projdatum->ey)*cos(projdatum->ex)+cos(projdatum->ez)*sin(projdatum->ex)); /*error here too many * */
    projdatum->r33=projdatum->m*cos(projdatum->ey)*cos(projdatum->ex);
    break;
  }
}
/* procedure to convert lat lon h from WGS84 to defined projection
*  The data is stored in the database in lat, lon , h WGS84 and need only to be converted in
*    FMaps->projdatum
*/
void gtk_fmaps_project(GtkFMaps *fmaps, gint tablenb, gdouble *lon, gdouble *lat, gdouble *h)
{
  gdouble x,y,z;

  g_return_if_fail(fmaps!=NULL);

  if (fmaps->projdatum->datumtransform!=DAT_NONE)
  {
    /* convert first in cartesian */
    gtk_fmaps_lonlath_2_xyz(fmaps->projdatum,(*lon*DEG_TO_RAD),(*lat*DEG_TO_RAD),*h,&x,&y,&z);

    /* convert in new datum*/
    switch(fmaps->projdatum->datumtransform)
    {
    case DAT_NONE:
        break;
    case DAT_MOLENDENSKY:
        gtk_fmaps_molendensky(fmaps->projdatum,TRUE,&x,&y,&z);
        break;
    case DAT_BURSA_WOLFE:
        gtk_fmaps_bursa_wolfe(fmaps->projdatum,TRUE,&x,&y,&z);
    }

    /* bring back in lon lat h */
    gtk_fmaps_xyz_2_lonlath(fmaps->projdatum,x,y,z,lon,lat,h);
    *lon=*lon*RAD_TO_DEG;
    *lat=*lat*DEG_TO_RAD;
  }
}

---