Sunday, July 22, 2012

ArcGis JavaScript API and Google Places API


1.     Introduction

Google exposes its database of places for use with GIS  applications. Places can be queried through a spatial query or through attribute queries. In this document I will only illustrate the use of the spatial query in Google places API.  Before you can use the API you must create a Google account and activate the Google API places API for your account.
The Google API used here is the one illustrated below.

More details about parameters can be found at the URL:

In this example I will only use two parameters, the location and the distance. Information returned is a list of maximum 60 places in trunks of 20 items. All results are returned in WGS 84 coordinates.

2.     Using Place API on a map with a base layer of 102100 SR

Using the spatial reference 102100 makes it very easy to convert the WGS 84 coordinates into map points of the ArcGis map. The result of this example looks like this on our map :

Integrating the places with an existing ArcGis JavaScript application can be done using the following steps:

1.       Start a point selection on the map returning a point.

2.       With the point , start a first call to the Google API with the location specified in WGS84 coordinates.

3.       Use the resulting array of maximum 20 results to add the contents to the graphic layer of the map. You can use the Google icons as symbol, and use the Google attributes for displaying of an info window

4.       If a next token is returned, you can use it to request  another 20 results. Repeat this until the next token is not defined. Current a maximum of 60 locations are returned.

In this implementation I used the XML version because I have already a cross domain AJAX method available for XML output. The Google places API implementation into the JavaScript framework looks like this:

findPlaces: function (point, radiusAPI, callback) {
  try {
     findPlacesCallback = callback;
     fullResults = new Array();
     location = point.y.toString() + "," + point.x.toString();
     radius = radiusAPI;
     queryPlacesByToken(null);
 } catch (e) {
     logging.logMessage("E", "error in Google place API interface ->" + err.message, "GisExtra/findPlaces");
 }
},

findPlacesByToken: function (token, callback) {
  try {
     findPlacesCallback = callback;
     queryPlacesByToken(token);
  } catch (e) {
     logging.logMessage("E", "error in Google place token API interface ->" + err.message, "GisExtra/findPlacesByToken");
  }
}

The first method is used to start the Google places search, the second is called when a next token is available in the result indicating that more results can be retrieved.  The actual initiation of the  Google API is done in the method ‘queryPlacesByToken’. A callback method must be provided by the calling method for returning the Google places data.

function queryPlacesByToken(token) {
   var data;
   if (token == null) {
     data = {
        location: location,
        sensor: false,
        radius: radius,
        rankBy: 'distance',
        key: GoogleAPIKey
     };
   } else {
     data = {
        sensor: false,
        pagetoken: token,
        key: GoogleAPIKey
     };
   }
   var url = googleUrlPlaces + "?" + $.param(data);
   var doneOnce = false;
   requestCrossDomain(url, function (result) {
     try {
       if (!doneOnce && findPlacesCallback != null) {
         doneOnce = true;
         var jsonContents = $.xml2json(result.results[0]);
         if (jsonContents.status == "OK") {
           findPlacesSuccess(jsonContents.result, jsonContents.next_page_token);
         } else {
           findPlacesSuccess(null, jsonContents.status);
         }
       }
     } catch (e) {
        logging.logMessage("E", "error in requestCrossDomain ->" + e.message, "GisExtra/queryPlacesByToken");
     }
   });
 }

When the result is returned from the Google API GET the status must be OK otherwise an error occurred. Beside the list of places, a next token is also returned to the calling method (function).
 I noticed that requesting to quickly for the next array of results from the Google places API, the same results were returned.  To solve this issue, I first let the calling method build the graphics for the map and only at the end call for the next array of results. A wait is also introduced until all graphics are add to the graphic layer before  a next array of places is retrieved. That is the reason why a token is returned and through the second method  ‘findPlacesByToken’ in the GIS library, the application is able to request the next list of places.

timeout = setTimeout(nextResult, 2000);

function nextResult() {
  if (GisOperation.getMap().graphics.graphics.length == googlePlacesCount) {
    if (nextToken) {
       GisExtra.findPlacesByToken(nextToken, processResults);
    } else {
    }
    clearTimeout(timeout);
  }
}

3.     Using Place API on a map with a base layer of EPSG:31370


The real problems comes when you use a projection based on another datum than WGS84.
After the selection of a point in step 1, you need to do a project transformation from 31370 to WGS84 because Google requires a point in latitude and longitude coordinates for doing the spatial query.

As from ArcGis 10.1 you can  add a transformation to the rest call to the project method of the geometry service.  In case of 31370 spatial reference you can use Belge_1972_To_WGS_1984_3 datum transformation. And this works fine.
However when using the inverse functionality of the datum transformation, the method gives me two problems.

First, when requesting the transformation of 20 points the rest service gives an error, using i.e. 15 points does returns results. This could be a proxy problem. This can be solved by a workaround through splitting the datum transformation into maximum 15 points, using maximum two calls to the rest service for transforming the points.
The biggest problem is that the results have a serious offset from the actual point locations.

So I started by creating a datum transformation in JavaScript, you can find more of this on the next sections. These are the results of one point.
Green point: correct location, based on a map projection of 102100.
Yellow point, map projection 31370 :     horizontal: 197m, vertical: 120m, distance: 214m
Bleu point, map projection 31370:          horizontal: 71m, vertical: 22m, distance: 74,5m

The yellow point is the result of the use of the method project of the geometry rest service. The bleu point corresponds to the custom JavaScript method used.  Still not an exciting result, but already a lot better by using the ESRI transformation. I put a thread on the ESRI JavaScript forum, but as for now no respond.


4.     Creating a project transformation WGS84 to EPSG:31370


The illustration below shows how we can create a custom projection transformation in JavaScript, client side.

The following steps has to be done for doing the project transformation from WGS84 to Lambert 72 (31370).

1.       Transform  geographic coordinates of the geodetic datum into geocentric coordinates.

2.       Perform a translation on the geocentric coordinates using  the Helmert 7-parameter transformations (position vector or coordinate frame).  There are other transformations possible depending on the required country projections.

3.       Transform geocentric coordinates back into geographic coordinates.

4.       Transform the geographic coordinates into the projection coordinates, in this example this is Lambert 1972 (31370).

There is a lot of literature that is useful for this kind of transformation. At the end you can find URL’s to some documents. A lot of these are purely mathematic formula’s, not exiting stuff for a developer, but lucky I studied mathematics, and this was a great help.

4.1 Doing transformations from geographic to geocentric coordinates

This kind of transformation depends on the eclipse (geodetic datum) used in the geometric coordinates and requires some parameter values. The nature of the required parameter looks like:

    // ellipse parameters
    CoordTransform.ellipse = {
      WGS84: { a: 6378137, b: 6356752.3142, f: 1 / 298.257223563, e: 0.081819191, e2: 0.0066943800158945 },
      Intl1924: { a: 6378388.000, b: 6356911.946, f: 1 / 297.0, e: 0.08199189, e2: 0.00672267002 }
    };

If  the geodetic datum is different of  the ones above, you will certainly find it on the internet.
Given a point with latitude and longitude, this is how you can create the geocentric coordinates.  Input is an eclipse e1. More details of this can be found on the URL’s at the end of the document.

        var lat = point.lat * Math.PI / 180;
         var lon = point.lng * Math.PI / 180;
        var a = e1.a, b = e1.b;
        var sinPhi = Math.sin(lat);
        var cosPhi = Math.cos(lat);
        var sinLambda = Math.sin(lon);
        var cosLambda = Math.cos(lon);
        var H = 200;  // standard hight
        var eSq = (a * a - b * b) / (a * a);
        var nu = a / Math.sqrt(1 - eSq * sinPhi * sinPhi);
        var x1 = (nu + H) * cosPhi * cosLambda;
        var y1 = (nu + H) * cosPhi * sinLambda;
        var z1 = ((1 - eSq) * nu + H) * sinPhi;

4.2 Translate the geocentric coordinates from one datum to another

For doing this I use the Helmert 7 parameter transformation. At the site of ESRI you can find a list of all defined transformations with the 7 parameters.

More information is available on the internet. 
Input are the 7 parameters of the transformation.  For Lambert 72 this is :

            // helmert transform parameters from WGS84 to other datums
            CoordTransform.datumTransform = {
                // Helmert parameters
                toLambert72: { tx: 106.868628, ty: -52.297783, tz: 103.723893,  //
                    rx: 0.33657, ry: -0.456955, rz: 1.842183, // sec
                    s: -1.2747,
               // Other parameters for step 4
                    latF: 1.57079633,
                    lngF: 0.07604294,
                    lat1Par: 0.86975574,
                    lat2Par: 0.89302680,
                    eastF: 150000.01256,
                    northF: 5400088.4378
                }
            }

The result is returned in x2,y2,z2

        var tx = t.tx, ty = t.ty, tz = t.tz;
        var rx = (t.rx / 3600).toRad();  // normalise seconds to radians
        var ry = (t.ry / 3600).toRad();
        var rz = (t.rz / 3600).toRad();
        var s1 = (t.s / 1e6) + 1;          // normalise ppm to (s+1)

        // apply transform
        var x2 = tx + s1 * (x1  - y1 * rz + z1 * ry);
        var y2 = ty + s1 * (y1  + x1 * rz - z1 * rx);
        var z2 = tz + s1 * (z1  - x1 * ry + y1 * rx);

This is a simple matrix transformation.

4.3 Doing transformations from geocentric to geographic  coordinates


This is the inverse operation of 4.1. Input is here the eclipse parameters from to target datum. In our case this is the International 1924 which is the datum of Lambert 72.

The transformation is based on the Browring Algorithm, an takes as input the eclipse parameters of the target eclipse. The result of the calculation is latitude lat2 and longitude ln2. The height H is not important.
        a = e2.a, b = e2.b;
        var precision = 2 / a;  // results accurate to around 2 metres
        eSq = (a * a - b * b) / (a * a);
        var p = Math.sqrt(x2 * x2 + y2 * y2);
        var lat2 = Math.atan2(z2, p * (1 - eSq)), phiP = 2 * Math.PI;
        while (Math.abs(lat2 - phiP) > precision) {
            nu = a / Math.sqrt(1 - eSq * Math.sin(lat2) * Math.sin(lat2));
            phiP = lat2;
            lat2 = Math.atan2(z2 + eSq * nu * Math.sin(lat2), p);
        }
        var lng2 = Math.atan2(y2, x2);
        H = p / Math.cos(lat2) - nu;

4.4  Transform the geographic coordinates into the projection coordinates

The calculation needed to create the projection coordinates out of geographic coordinates depends on the projection. What is illustrated here is only valid for the projection Lambert 72 (31370) but can you give some indication how to do the transformation.
As mentioned earlier, some extra parameter values has been added to the transformation parameters of Lambert 72.
Transformation from geographic coordinates into projection coordinates are based on a method that is specified in the same table as you could find the Helmert parameters (see 4.2). In our case you will see the value CF which means coordinate frame. In our case we will use. The full name of the transformation is regular Lambert Conic Conformal (2SP Belgium) .
Input is a point with latitude and longitude, eclipse and transformation parameters. The transformation parameters are those added after the Helmert parameters. The result are the north and east values of a point in Lambert 72.

    function conicConformalSP2(pointLatLng, e, t) {
        var m1 = Math.cos(t.lat1Par) / Math.sqrt((1 - e.e2 * Math.sin(t.lat1Par) * Math.sin(t.lat1Par)));
        var m2 = Math.cos(t.lat2Par) / Math.sqrt((1 - e.e2 * Math.sin(t.lat2Par) * Math.sin(t.lat2Par)));
        var t1 = Math.tan(Math.PI / 4 - t.lat1Par / 2) / Math.pow((1 - e.e * Math.sin(t.lat1Par)) / (1 + e.e * Math.sin(t.lat1Par)), e.e / 2);
        var t2 = Math.tan(Math.PI / 4 - t.lat2Par / 2) / Math.pow((1 - e.e * Math.sin(t.lat2Par)) / (1 + e.e * Math.sin(t.lat2Par)), e.e / 2);
        var tF = Math.tan(Math.PI / 4 - t.latF / 2) / Math.pow((1 - e.e * Math.sin(t.latF)) / (1 + e.e * Math.sin(t.latF)), e.e / 2);
        if (Math.abs(tF) < 0.0000001)
            tF = 0;
        var tPt = Math.tan(Math.PI / 4 - pointLatLng.lat / 2) / Math.pow((1 - e.e * Math.sin(pointLatLng.lat)) / (1 + e.e * Math.sin(pointLatLng.lat)), e.e / 2);
        var n = (Math.log(m1) - Math.log(m2)) / (Math.log(t1) - Math.log(t2));
        var F = m1 / (n * Math.pow(t1, n));
        var r = e.a * F * Math.pow(tPt, n);
        var rF = 0;
        if (tF > 0)
            rF = e.a * F * Math.pow(tF, n);
        var tetra = n * (pointLatLng.lng - t.lngF);
        var a = 0.00014204;
        var x = t.eastF + r * Math.sin(tetra - a);
        var y = t.northF + rF - r * Math.cos(tetra - a);
        return { x: x, y: y };
    }

5.     Conclusion

As you saw earlier, doing the transformation in code instead of doing  a project with the geometry services results in better accuracy but still is a large difference with the actual position. Only using a map with as spatial reference 102100 results in a correct location. I don’t know why there is still a visible offset of the real location, I suspect that the Helmert transformation lacks precision. There exist also a 10 parameter transformation ‘Molodensky -Badekas method’ but I could not find the parameters for Lambert 72.
The use of JavaScript transformation with greater accuracy solves a lot of asynchronous handling of data that is needed in case geometry services are used.
I cannot explain why from Lambert 72 to WGS84 the transformation in the geometry service is correct and in the other direction gives a large offset.

6.     Interesting links

Example of datum conversion for WGS84 to OSGB36 (Great-Britain).  Some JavaScript code is based on the example.


Theoretical explanation the different eclipses that is used for the geographic coordinates and different transformation towards projections.


Visual explanation about geographic coordinates , translation and formulas.


Eclipse parameters for different datums.



Wednesday, July 4, 2012

Using different projections within ArcGis JavaScript Applications


1.     Introduction


Did you ever have a displacement of points of 100m on your map, then you clearly run into problems with datum corrections. This document explains how I tackle this issue , and how ArcGis 10.1 helps you.

Public geo databases as ‘Openstreet map’  , ‘Google maps’  or ‘Google Streetview’ are all based on the spatial reference WGS84 and this is in most cases not the spatial reference of the map.  When trying to integrate this information with a custom ArcGis JavaScript, you will need to project  your data towards WGS 84 spatial reference coordinates.

To project your data into WGS 84 spatial reference, ESRI expose a project transformation service through the method project of the  geometry service of an ArcGis server.  As from ArcGis 10.1 and JavaScript API 3.0 the project method has added a transformation parameter. The list of all transformations can be found at the URL  http://help.arcgis.com/en/webapi/javascript/arcgis/help/jshelp/dattrans.htm

In pre version 3.0, if a transformation is required, a SOAP  web service must be used, however using this in JavaScript requires some extra data processing,  a JSON approach is more practical.

Sometimes you will need using more than one transformation in the case no direct transformation of projection to another projection is available as outlined in the figure below.


2.     Transformation between WGS84 and WebMercator (102100/3857)


Because I could not find a transformation in the list on the URL specified in the ‘Introduction’ I wrote two  functions that handle this transformations in both directions based on some C++ code I found. Doing so, I did not need an async REST call to the ArcGIs Geometry service, speeding up the transformation.

The code for doing these transformations are :

function toGeographic(xMercator, yMercator) {
  if (Math.abs(xMercator) < 180 && Math.abs(yMercator) < 90)
     return null;
  if ((Math.abs(xMercator) > 20037508.3427892) || (Math.abs(yMercator) > 20037508.3427892))
      return null;
  var x = xMercator;
  var y = yMercator;
  var w1 = x = x / 6378137.0;
  var w2 = x * 57.295779513082323;
  var w3 = Math.floor((x + 180.0) / 360.0);
  x = w2 - (w3 * 360.0);
  y = (1.5707963267948966 - (2.0 * Math.atan(Math.exp((-1.0 * y) / 6378137.0)))) * 57.295779513082323;
  return {
    x: x,
    y: y
  }
}

function toWebMercator(xlon, ylat) {
  if ((Math.abs(mercatorX_lon) > 180 || Math.abs(mercatorY_lat) > 90))
    return null;
  var num = xlon * 0.017453292519943295;
  var x = 6378137.0 * num;
  var a = ylat * 0.017453292519943295;
  var y = 3189068.5 * Math.Log((1.0 + Math.Sin(a)) / (1.0 - Math.Sin(a)));
  return {
     x: x,
     y: y
  }
}

3.     Other Project transformation


Creation of a generic  function for doing projection transformation requires as input either a wkid or wkt of the transformation that corresponds to an item of the transformation list found at the URL specified earlier.  In a real world application a transformation happens between an input spatial reference (wkid) and an output spatial reference (wkid).
I could not find a relation between the transformation list wkid’s and the wkid’s  of spatial references.
The only solution for this lack of automatic solution based on the two wkid’s was to add a parameter at the layer configuration specifying the transformation needed for the WGS 84 spatial references conversion.  This solves all my issues for the transformation of country spatial references to WGS 84 spatial reference, making integration with Google Streetview possible and other WGS 84 related geographical data.
I also had a spatial reference problem in case WMS layers are used for doing spatial queries. A WMS query requires the same spatial reference as the WMS service expose. A solution is to add a parameter to the WMS configuration specifying the transformation needed for modifying the base layer spatial reference to the WMS spatial reference for the geometry involved.
The result of the project conversion consist of two cases, one that looks if it involves a Webmercator / WGS 84 transformation, and a second that will use the geometry service of the ArcGis server. The method requires version 3.0 of the ArcGis JavaScript API and ArcGis Server 10.1.

convertProjection: function (points, toSpacialReference, transformation, callback) {
   var wkidIn = points[0].spatialReference.wkid;
   var wkidOut = toSpacialReference.wkid;
   var result = null;
   var resultPoints = new Array();
   var point;
   if ((wkidIn == '102100' || wkidIn == '3857') && wkidOut == '4326') {
     for (var i = 0; i < points.length; i++) {
        result = toGeographic(points[0].x, points[0].y);
        point = new esri.geometry.Point([result.x, result.y], new esri.SpatialReference({ wkid: 4326 }));
        resultPoints.push(point);
     }
     callback(resultPoints);
   }
   else if (wkidIn == '4326' && (wkidOut == '102100' || wkidIn == '3857')) {
     for (var i = 0; i < points.length; i++) {
        result = toWebMercator(points[0].x, points[0].y);
        point = new esri.geometry.Point([result.x, result.y], new esri.SpatialReference({ wkid: 102100 }));
        resultPoints.push(point);
     }
     callback(resultPoints);
   }
   else {
     var params = new esri.tasks.ProjectParameters();
     params.geometries = points;
     params.outSR = toSpacialReference;
     params.transformation = transformation;
     params.transformationForward = false;
     GisOperation.getGeometryService().project(params,
       function (projectedPoints) {
          callback(projectedPoints);
       }, function errorprojected(err) {
          logging.logMessage("E", "projection convertion failed -->" + err.message, "GisGeoProcessing");
          callBack(null);
     });
   }
}

4.     Transformation in an application


The above transformation worked well for me having as spatial reference of Lambert 72 (31370) of the map, and using Google Streetview to show the location selected. The transformation used was 1610.