Commit 875a8689 authored by Richard Didier's avatar Richard Didier
Browse files

Support for nadgrids : It is a port of the CTABLE structure from PROJ.4 in the... 

Support for nadgrids : It is a port of the CTABLE structure from PROJ.4 in the trunk. See ticket #81
parent 6e1ba703

lib/defs/ntf_r93.gsb.js

0 → 100644
+1 −0

File added.

Preview size limit exceeded, changes collapsed.

lib/proj4js.js

+266 −6
Original line number Diff line number Diff line
@@ -153,13 +153,53 @@ var Proj4js = { 
          return point;

      }



      //DGR: 2012-07-29 : add nadgrids support (begin)

      var src_a = source.a;

      var src_es = source.es;



      var dst_a = dest.a;

      var dst_es = dest.es;



      var fallback= source.datum_type;

      // If this datum requires grid shifts, then apply it to geodetic coordinates.

      if( fallback == Proj4js.common.PJD_GRIDSHIFT )

      {

          if (this.apply_gridshift( source, 0, point )==0) {

            src_a = Proj4js.common.SRS_WGS84_SEMIMAJOR;

            src_es = Proj4js.common.SRS_WGS84_ESQUARED;

          } else {



              // try 3 or 7 params transformation or nothing ?

              if (!source.datum_params) {

                  return point;

              }

              var wp= 1.0;

              for (var i= 0, l= source.datum_params.length; i<l; i++) {

                wp*= source.datum_params[i];

              }

              if (wp==0.0) {

                return point;

              }

              fallback= source.datum_params.length>3?

                Proj4js.common.PJD_7PARAM

              : Proj4js.common.PJD_3PARAM;

              // CHECK_RETURN;

          }

      }



      if( dest.datum_type == Proj4js.common.PJD_GRIDSHIFT )

      {

          dst_a = Proj4js.common.SRS_WGS84_SEMIMAJOR;

          dst_es = Proj4js.common.SRS_WGS84_ESQUARED;

      }

      // Do we need to go through geocentric coordinates?

      if( source.es != dest.es || source.a != dest.a

          || source.datum_type == Proj4js.common.PJD_3PARAM

          || source.datum_type == Proj4js.common.PJD_7PARAM

      if( src_es != dst_es || src_a != dst_a

          || fallback == Proj4js.common.PJD_3PARAM

          || fallback == Proj4js.common.PJD_7PARAM

          || dest.datum_type == Proj4js.common.PJD_3PARAM

          || dest.datum_type == Proj4js.common.PJD_7PARAM)

      {

      //DGR: 2012-07-29 : add nadgrids support (end)



        // Convert to geocentric coordinates.

        source.geodetic_to_geocentric( point );
@@ -181,9 +221,70 @@ var Proj4js = { 
          // CHECK_RETURN;

      }



      // Apply grid shift to destination if required

      if( dest.datum_type == Proj4js.common.PJD_GRIDSHIFT )

      {

          this.apply_gridshift( dest, 1, point);

          // CHECK_RETURN;

      }



      return point;

    }, // cs_datum_transform



    /**

     * This is the real workhorse, given a gridlist

     * DGR: 2012-07-29 addition based on proj4 trunk

     */

    apply_gridshift : function(srs,inverse,point) {

        if (srs.grids==null || srs.grids.length==0) {

            return -38;

        }

        var input= {"x":point.x, "y":point.y};

        var output= {"x":Number.NaN, "y":Number.NaN};

        /* keep trying till we find a table that works */

        var onlyMandatoryGrids= false;

        for (var i= 0, l= srs.grids.length; i<l; i++) {

            var gi= srs.grids[i];

            onlyMandatoryGrids= gi.mandatory;

            var ct= gi.grid;

            if (ct==null) {

                if (gi.mandatory) {

                    this.reportError("unable to find '"+gi.name+"' grid.");

                    return -48;

                }

                continue;//optional grid

            } 

            /* skip tables that don't match our point at all.  */

            var epsilon= (Math.abs(ct.del[1])+Math.abs(ct.del[0]))/10000.0;

            if( ct.ll[1]-epsilon>input.y || ct.ll[0]-epsilon>input.x ||

                ct.ll[1]+(ct.lim[1]-1)*ct.del[1]+epsilon<input.y ||

                ct.ll[0]+(ct.lim[0]-1)*ct.del[0]+epsilon<input.x ) {

                continue;

            }

            /* If we have child nodes, check to see if any of them apply. */

            /* TODO : only plain grid has implemented ... */

            /* we found a more refined child node to use */

            /* load the grid shift info if we don't have it. */

            /* TODO : Proj4js.grids pre-loaded (as they can be huge ...) */

            output= Proj4js.common.nad_cvt(input, inverse, ct);

            if (!isNaN(output.x)) {

                break;

            }

        }

        if (isNaN(output.x)) {

            if (!onlyMandatoryGrids) {

                this.reportError("failed to find a grid shift table for location '"+

                    input.x*Proj4js.common.R2D+" "+input.y*Proj4js.common.R2D+

                    " tried: '"+srs.nadgrids+"'");

                return -48;

            }

            return -1;//FIXME: no shift applied ...

        }

        point.x= output.x;

        point.y= output.y;

        return 0;

    },



    /**

     * Function: adjust_axis

     * Normalize or de-normalized the x/y/z axes.  The normal form is "enu"
@@ -942,6 +1043,7 @@ Proj4js.Proj = Proj4js.Class({ 
                                this.axis= paramVal;

                             } //FIXME: be silent ?

                             break;

              case "wktext": break;//DGR 2012-07-29

              case "no_defs": break;

              default: //alert("Unrecognized parameter: " + paramName);

          } // switch()
@@ -957,7 +1059,31 @@ Proj4js.Proj = Proj4js.Class({ 
 */

  deriveConstants: function() {

      // DGR 2011-03-20 : nagrids -> nadgrids

      if (this.nadgrids == '@null') this.datumCode = 'none';

      if (this.nadgrids && this.nadgrids.length==0) {

          this.nadgrids= null;

      }

      if (this.nadgrids) {

        if (this.nadgrids == '@null' && !this.datum_params) {

            this.datumCode = 'none';

        } else {

          this.grids= this.nadgrids.split(",");

          var g= null;

          for (var i= 0, l= this.grids.length; i<l; i++) {

            g= this.grids[i];

            var fg= g.split("@");

            this.grids[i]= {

              mandatory: fg.length==1,//@=> optional grid (no error if not found)

              name:fg[fg.length-1],

              grid: Proj4js.grids[fg[fg.length-1]]//FIXME: grids loading ...

            };

            if (this.grids[i].mandatory && !this.grids[i].grid) {

              Proj4js.reportError("Missing '"+this.grids[i].name+"'");

            }

          }

        }

        // DGR, 2011-03-20: grids is an array of objects that hold

        // the loaded grids, its name and the mandatory informations of it.

      }

      if (this.datumCode && this.datumCode != 'none') {

        var datumDef = Proj4js.Datum[this.datumCode];

        if (datumDef) {
@@ -1058,6 +1184,7 @@ Proj4js.common = { 
  PJD_WGS84    : 4,   // WGS84 or equivalent

  PJD_NODATUM  : 5,   // WGS84 or equivalent

  SRS_WGS84_SEMIMAJOR : 6378137.0,  // only used in grid shift transforms

  SRS_WGS84_ESQUARED : 0.0066943799901413165, //DGR: 2012-07-29



  // ellipoid pj_set_ell.c

  SIXTH : .1666666666666666667, /* 1/6 */
@@ -1315,6 +1442,110 @@ Proj4js.common = { 
    return phi;

  },



  /**

   * Determine correction values

   * source: nad_intr.c (DGR: 2012-07-29)

   */

  nad_intr: function(pin,ct) {

    // force computation by decreasing by 1e-7 to be as closed as possible

    // from computation under C:C++ by leveraging rounding problems ...

    var t= {"x":(pin.x-1.e-7)/ct.del[0],"y":(pin.y-1e-7)/ct.del[1]};

    var indx= {"x":Math.floor(t.x),"y":Math.floor(t.y)};

    var frct= {"x":t.x-1.0*indx.x,"y":t.y-1.0*indx.y};

    var val= {"x":Number.NaN,"y":Number.NaN};

    var inx;

    if (indx.x<0) {

      if (!(indx.x==-1 && frct.x>0.99999999999)) {

        return val;

      }

      ++indx.x;

      frct.x= 0.0;

    } else {

      inx= indx.x+1;

      if (inx>=ct.lim[0]) {

        if (!(inx==ct.lim[0] && frct.x<1e-11)) {

          return val;

        }

        --indx.x;

        frct.x= 1.0;

      }

    }

    if (indx.y<0) {

      if (!(indx.y==-1 && frct.y>0.99999999999)) {

        return val;

      }

      ++indx.y;

      frct.y= 0.0;

    } else {

      inx= indx.y+1;

      if (inx>=ct.lim[1]) {

        if (!(inx==ct.lim[1] && frct.y<1e-11)) {

          return val;

        }

        --indx.y;

        frct.y= 1.0;

      }

    }

    inx= (indx.y*ct.lim[0])+indx.x;

    var f00= {"x":ct.cvs[inx][0], "y":ct.cvs[inx][1]};

    inx++;

    var f10= {"x":ct.cvs[inx][0], "y":ct.cvs[inx][1]};

    inx+= ct.lim[0];

    var f11= {"x":ct.cvs[inx][0], "y":ct.cvs[inx][1]};

    inx--;

    var f01= {"x":ct.cvs[inx][0], "y":ct.cvs[inx][1]};

    var m11= frct.x*frct.y,             m10= frct.x*(1.0-frct.y),

        m00= (1.0-frct.x)*(1.0-frct.y), m01= (1.0-frct.x)*frct.y;

    val.x= (m00*f00.x + m10*f10.x + m01*f01.x + m11*f11.x);

    val.y= (m00*f00.y + m10*f10.y + m01*f01.y + m11*f11.y);

    return val;

  },



  /**

   * Correct value

   * source: nad_cvt.c (DGR: 2012-07-29)

   */

  nad_cvt: function(pin,inverse,ct) {

    var val= {"x":Number.NaN, "y":Number.NaN};

    if (isNaN(pin.x)) { return val; }

    var tb= {"x":pin.x, "y":pin.y};

    tb.x-= ct.ll[0];

    tb.y-= ct.ll[1];

    tb.x= Proj4js.common.adjust_lon(tb.x - Proj4js.common.PI) + Proj4js.common.PI;

    var t= Proj4js.common.nad_intr(tb,ct);

    if (inverse) {

      if (isNaN(t.x)) {

        return val;

      }

      t.x= tb.x + t.x;

      t.y= tb.y - t.y;

      var i= 9, tol= 1e-12;

      var dif, del;

      do {

        del= Proj4js.common.nad_intr(t,ct);

        if (isNaN(del.x)) {

          this.reportError("Inverse grid shift iteration failed, presumably at grid edge.  Using first approximation.");

          break;

        }

        dif= {"x":t.x-del.x-tb.x, "y":t.y+del.y-tb.y};

        t.x-= dif.x;

        t.y-= dif.y;

      } while (i-- && Math.abs(dif.x)>tol && Math.abs(dif.y)>tol);

      if (i<0) {

        this.reportError("Inverse grid shift iterator failed to converge.");

        return val;

      }

      val.x= Proj4js.common.adjust_lon(t.x+ct.ll[0]);

      val.y= t.y+ct.ll[1];

    } else {

      if (!isNaN(t.x)) {

          val.x= pin.x - t.x;

          val.y= pin.y + t.y;

      }

    }

    return val;

  },



/* meridinal distance for ellipsoid and inverse

**	8th degree - accurate to < 1e-5 meters when used in conjuction

**		with typical major axis values.
@@ -1363,11 +1594,19 @@ Proj4js.datum = Proj4js.Class({ 
        }

      }

    }

    // DGR 2011-03-21 : nadgrids support

    this.datum_type = proj.grids?

      Proj4js.common.PJD_GRIDSHIFT

    : this.datum_type;



    this.a = proj.a;    //datum object also uses these values

    this.b = proj.b;

    this.es = proj.es;

    this.ep2 = proj.ep2;

    this.datum_params = proj.datum_params;

    if (this.datum_type==Proj4js.common.PJD_GRIDSHIFT) {

      this.grids= proj.grids;

    }

  },



  /****************************************************************/
@@ -1394,8 +1633,10 @@ Proj4js.datum = Proj4js.Class({ 
              && this.datum_params[6] == dest.datum_params[6]);

    } else if ( this.datum_type == Proj4js.common.PJD_GRIDSHIFT ||

                dest.datum_type == Proj4js.common.PJD_GRIDSHIFT ) {

      alert("ERROR: Grid shift transformations are not implemented.");

      return false

      //alert("ERROR: Grid shift transformations are not implemented.");

      //return false

      //DGR 2012-07-29 lazy ...

      return this.nadgrids == dest.nadgrids;

    } else {

      return true; // datums are equal

    }
@@ -1889,4 +2130,23 @@ Proj4js.wktProjections = { 
  "Universal Transverse Mercator System": "utm"

};



// Based on proj4 CTABLE  structure :

// FIXME: better to have array instead of object holding longitudes, latitudes members

//        In the former case, one has to document index 0 is longitude and

//        1 is latitude ...

//        In the later case, grid object gets bigger !!!!

//        Solution 1 is chosen based on pj_gridinfo.c

Proj4js.grids= {

    "null":{                                // name of grid's file

        "ll": [-3.14159265, -1.57079633],   // lower-left coordinates in radians (longitude, latitude):

        "del":[ 3.14159265,  1.57079633],   // cell's size in radians (longitude, latitude):

        "lim":[ 3, 3],                      // number of nodes in longitude, latitude (including edges):

        "count":9,                          // total number of nodes

        "cvs":[                             // shifts : in ntv2 reverse order : lon, lat in radians ...

            [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], // for (lon= 0; lon<lim[0]; lon++) {

            [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], //   for (lat= 0; lat<lim[1]; lat++) { p= cvs[lat*lim[0]+lon]; }

            [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]  // }

        ]

    }

};