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.



1 comment:

  1. John thanks for this informative post. I have been banging my head against the wall after encountering the same offset problem with points coming out of NAD83 State Plane system. I look forward to reading the linked resources and this post in more detail.

    ~Seth

    ReplyDelete