Distance calculations using Postgis

Question

I have a distance calculation in Postgres on a site that I did not build.

I programmed the script that inserts new locations in php because I know php better the ruby.

I am trying to figure out how I can insert new geometry points for locations.

I have tried using the following but I get like 8000+ miles between the points...

The current points in the database works with distance... but any new points I try to insert get the wrong distances. Even if I add 2 new points to check:

I am using the following code with a csv file of city state and zipcodes the rebuild the database:

SAMPLE File:

2832    Hope Valley RI
2833    Hopkinton   RI
2835    Jamestown   RI
2836    Kenyon  RI
2837    Little Compton  RI
2838    Manville    RI
2839    Mapleville  RI
2840    Newport RI
2841    Newport RI
2842    Middletown  RI
2852    North Kingstown RI
2854    North Kingstown RI
2857    North Scituate  RI
2858    Oakland RI
2859    Pascoag RI
2860    Pawtucket   RI
2861    Pawtucket   RI
2862    Pawtucket   RI
2863    Central Falls   RI
2864    Cumberland  RI
2865    Lincoln RI
2871    Portsmouth  RI
2872    Prudence Island RI
2873    Rockville   RI
2874    Saunderstown    RI
2875    Shannock    RI
2876    Slatersville    RI
2877    Slocum  RI
2878    Tiverton    RI
2879    Wakefield   RI
2880    Wakefield   RI
2881    Kingston    RI
2882    Narragansett    RI
2883    Peace Dale  RI
2885    Warren  RI
2886    Warwick RI
2887    Warwick RI

My PHP insert Code

$content = file($selectfile1);

$posted_content = array();
list($rownum, $row) = each($content);
$posted_content[0] = explode(",", $row);
array_push($posted_content[0], "ID");
$count = 0;
$error = 0;
// iterate each row (1 post)
while (list($rownum, $row) = each($content))
{
    $count++;
    // extract fields from row columns
    $items = explode(",", $row);
    list($zip, $city1, $state1) = $items;       
    array_push($posted_content, $items);

    $getHTML ="http://dev.virtualearth.net/REST/v1/Locations/" . $zip . "?o=xml&key=AvtYnvs3UjKaIPdG5v1YaVBL_5-Rhg_zgUwoQgvTiTS9dMxJSreIanWVLvTzQc86";
    $getHTML1 ="http://dev.virtualearth.net/REST/v1/Locations/" . $city . "," .$state . "?o=xml&key=AvtYnvs3UjKaIPdG5v1YaVBL_5-Rhg_zgUwoQgvTiTS9dMxJSreIanWVLvTzQc86";
    // get the response from the Locations API and store it in a string
    $output = file_get_contents($getHTML);

    // create an XML element based on the XML string  
    $response = new SimpleXMLElement($output);

    // Extract data (e.g. latitude and longitude) from the results
    $lat = $response->ResourceSets->ResourceSet->Resources->Location->Point->Latitude;
    $lng = $response->ResourceSets->ResourceSet->Resources->Location->Point->Longitude;

    $lat1 = RTRIM($lat);
    $lat2 = LTRIM($lat1);
    $lng1 = RTRIM($lng);
    $lng2 = LTRIM($lng1);

    //$getJSON = "http://maps.googleapis.com/maps/api/geocode/json?address=" . $zip . "&sensor=false";
    //$contentJSON = file_get_contents($getJSON);
    //$Geocode_array = json_decode($contentJSON, true);

    //$lat = $Geocode_array['results'][0]['geometry']['location']['lat'];
    //$lng = $Geocode_array['results'][0]['geometry']['location']['lng'];

    reset($response);
    if(empty($lat2) || empty($lng2))
    {
        echo "<b>Error:</b>" . $city1 . ", " . $state1 . " " . $zip . "<br>";
        $error++;
        continue;
    }
    $vals = "ST_GeographyFromText('POINT(" . $lat2 . " " . $lng2 . ")')";

    $sql = "UPDATE locations SET coords =". $vals . " where zip ='" . $zip . "'";
    $rs = $conn->Execute($sql);
    if ($rs === false) die('Insert error: ' . $city . ', ' . $state . ' ' . $conn->ErrorMsg() . " SQL: " . $sql);


}
echo $error;
reset($content);
unset($GLOBALS["_SESSION"]["szContent"]);

I have tried the following and none have worked right with distance calculations (calculation code is below in the ruby coding)

$vals = "ST_GeographyFromText('POINT(" . $lat2 . " " . $lng2 . "), 4326')";


$vals = "ST_GeographyFromText('SRID=4326; POINT(" . $lat2 . " " . $lng2 . ")')";

Here is my database table that it is being inserted in:

CREATE TABLE locations
(
  id serial NOT NULL,
  coords geography(Point,4326),
  city character varying(255),
  state character varying(255),
  zip character varying(255),
  created_at timestamp without time zone NOT NULL,
  updated_at timestamp without time zone NOT NULL,
  cs character varying(255),
  alt text,
  CONSTRAINT locations_pkey PRIMARY KEY (id)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE locations
  OWNER TO postgres;
GRANT ALL ON TABLE locations TO postgres;
GRANT ALL ON TABLE locations TO loadmax;

-- Index: index_locations_on_coords

-- DROP INDEX index_locations_on_coords;

CREATE INDEX index_locations_on_coords
  ON locations
  USING gist
  (coords);

Ruby coding that calculates the distance on my ruby site...

class Location < ActiveRecord::Base
  attr_accessible :city, :coords, :state, :zip
  set_rgeo_factory_for_column(:coords,  RGeo::Geographic.spherical_factory(:srid => 4326))

  def self.validate_cs_and_return_proper(location)
    return "ANYWHERE" if location.upcase == "ANYWHERE" || location.nil? || location == ""
    state = location[-2,2]
    city = location.gsub(%r(, [a-zA-Z]+), '').strip 
    l = first(:conditions => {:city => city.downcase, :state => state.downcase}) 
    if l.nil?
      new_l = where("(levenshtein(city, '#{city}') <= 4) AND state = '#{state}'").order("levenshtein(city, '#{city}') ASC").first
      return "#{new_l.city}, #{new_l.state}"
    else
      return location
    end
  end

  def self.city_state_exists?(location)
    return true if location.upcase == "ANYWHERE" || location.nil? || location == ""
    state = location[-2,2]
    city = location.gsub(%r(, [a-zA-Z]+), '').strip 
    l = first(:conditions => {:city => city.downcase, :state => state.downcase}) 
    if l.nil?
      new_l = where("(levenshtein(city, '#{city}') <= 4) AND state = '#{state}'").order("levenshtein(city, '#{city}') ASC").first
      return true unless new_l.blank?
      return false
    else
      return true
    end
  end

  def self.to_point(location)
    return location unless location.class == String
    return nil if location.upcase=='ANYWHERE'
    state = location[-2,2]
    city = location.gsub(%r(, [a-zA-Z]+), '').strip 
    l = first(:conditions => {:city => city.downcase, :state => state.downcase}) 
    l.coords
  end

  def self.to_cs(point)
    return "ANYWHERE" if point.nil? || (point.x == 0.0 && point.y == 0.0)
    l = first(:conditions => {:coords => point})
    l[:state] = l[:state].upcase
    c_parts = l[:city].split(" ")
    c_parts.each {|part| part.capitalize! }
    l[:city] = c_parts.join(" ") 
    return "#{l[:city]}, #{l[:state].upcase}"
  end

  def self.city_state(point)
    return "ANYWHERE" if point.nil? || (point.x == 0.0 && point.y == 0.0)
    city_state = {}
    l = first(:conditions => {:coords => point})
    city_state[:city] = l[:city]
    city_state[:state] = l[:state]
    city_state
  end

  def self.distance_between(locale_one, locale_two)
    #locale_one and locale_two must be points
    select = "ST_Distance(gg1, gg2)"
    from = "(SELECT ST_GeographyFromText('#{locale_one}') AS gg1, ST_GeographyFromText('#{locale_two}') AS gg2) AS foo"
    distance = Location.select(select).from(from)
    #1609.344 is how many meters are in a mile
    "%.2f" % (distance[0].st_distance.to_f/1609.344)
  end

  def self.check_unknown(cs)
    report = Report.first( :conditions => { :metric => "unknown_locations", :name => cs } )
    #i.e. if it's been tried 25 times, add it to the locations table.
    if report.value > 5
      city, state = cs.split(",")
      state = state.strip
      google_api_results = HTTParty.get("http://maps.googleapis.com/maps/api/geocode/json?address=#{cs.gsub(/ /, "%20")}&sensor=true")
      geo_data = google_api_results.parsed_response
      lat = geo_data["results"][0]["geometry"]["location"]["lat"]
      lng =  geo_data["results"][0]["geometry"]["location"]["lng"]
      l = Location.new({ 
        :city => city,
        :state => state,
        :coords => "ST_GeographyFromText('POINT( #{lat} #{lng} )')", 
        :cs => "#{city}, #{state}"
      })
      l.save
      Report.delete(report.id)
      return l
    else
      return false
    end
  end

end

Answer 1

The axis order is important. If one assumes the input is (lat lng), (and I don't blame most people for thinking that way,) you will get unexpected results:

SELECT ST_AsLatLonText('POINT(37.215396881103516 -93.295166015625)');

86°42'17.402"S 142°47'4.571"W or La Gorce Mountains, Antarctica

And swapping them around to (lng lat):

SELECT ST_AsLatLonText('POINT(-93.295166015625 37.215396881103516)');

37°12'55.429"N 93°17'42.598"W or 400-504 W Chestnut Expy Springfield, MO 65801, USA

Use the correct axis order for PostGIS, which is (X Y) or (long lat). Therefore use:

$vals = "ST_GeographyFromText('POINT(" . $lng2 . " " . $lat2 . ")')"