Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

So I have created a polygon using Google Earth pro. I've exported it into a kml file, which I can then parse using ruby's XML/HTTP-parsing Nokogiri gem into an array of WKT points (LinearRing) thats soon to be a polygon. I'm feeding the geographic points into RGeo to try to get the equivalent points to make up a projected polygon originating out of the conversion from the geographic polygon.

class ReadKmlCoords
  attr_reader :ary, :coords
def initialize
    proj4 = "+proj=lcc +lat_1=34.83333333333334 +lat_2=32.5 +lat_0=31.83333333333333 +lon_0=-81 +x_0=609600 +y_0=0 +ellps=GRS80 +to_meter=0.3048 +no_defs"
    @factory = RGeo::Geographic.projected_factory( :projection_proj4 => proj4, :projection_srid => 3361, :buffer_resolution => 8 )
    @flat = RGeo::CoordSys::Proj4.create(proj4)
    @curved = @flat.get_geographic
  end

  def go(fn=nil)
    fn = fn || "gvilleblob.kml"
    x = Nokogiri::XML(File.open(fn).read)
    @ary = x.css("//coordinates").children.first.inner_text.gsub("\n", "").gsub("\t", "").split(",0")
    @ary.pop #removes spurious last datapoint
    @coords = []
    r = RGeo::WKRep::WKTParser.new
    @ary.each do |pair|
      lon = pair.split(",").first
      lat = pair.split(",").last
      geographic_point = r.parse("POINT (#{lon} #{lat})")
      @coords << RGeo::CoordSys::Proj4.transform(@curved, geographic_point, @flat, @factory)
    end
    @factory.polygon(@factory.projection_factory.linear_ring(@coords))
  end
end

This polygon is then to be stored in my database, which, in its projected form, will allow me to do PostGIS functions like st_overlaps, thus finding all the parcels lying under my big blob just made in Google Earth. You see, it seems there is some limitation on doing 'spatial joins' when using geographic topologies...right now I'm only able to get the which_parcels_are_under_this_manmade_blob method (PostGIS's st_overlaps method, RGeo's #overlaps? method, etc) to work when using the flat, projected coordinates. PostGIS scans about 100 parcels a second against another polygon vs RGeo's 1 parcel/s. So whatever I'm doing, it has to work with PostGIS.

Anyway, back to the processing, Unfortunately, the polygon produced by this code is inaccurate.

#An example of the (inaccurate, supposed-to-be-projected) polygon's points
#<RGeo::Geos::CAPIPolygonImpl:0x3fc5e22774a4 "POLYGON ((1577662.192624186 85.0511287, 1577366.1694666024 85.0511287, 1576894.8961126604 85.0511287, 1576139.5482155932 85.0511287, 1575483.7687068093 85.0511287

#And here are the first couple geographic coordinates (pre-WKT), the `@ary` in the code above 
@ary = ["-82.40799896121337,34.87381058338332", " -82.40898466663725,34.87373538049771", " -82.41054451142325,34.87304749870849", " -82.41304553437641,34.87200437474734", " -82.41521498428749,34.87098775443968", " -82.41738128573057,34.86997258001951", ...] 

I've used QGIS to check out the datapoints and it turns out the left number of each pair seems to be the correct projected x axis coordinate, while the right number is inaccurate.

ALTERNATIVELY, as a workaround solution, if there is any way I could do a an overlaps?/st_overlaps equivalent function using PostGIS which works with geographic coordinates, I'm all ears to that as well. Again, my purpose isn't to convert geographic coordinates in projected coordinates. I simply want to determine which parcels fall inside of a polygon, which originally is drawn in Google Earth. I am able to convert the projected points into geographic points, just not vice versa.

share|improve this question
    
It would be useful to see an example of the input that produced the output you show. –  whuber Aug 15 '13 at 22:43
    
You got me, I said the latitude was being kept, but that's not right at all. I was way off. Everything else is true so far, the big numbers in the projected coordinate set appear to be accurate. –  boulder_ruby Aug 15 '13 at 23:00
    
Definitely not good Y values there. I suspect lat isn't being set correctly somewhere. I guess you've confirmed that you're setting a good lat value in your @ary.each do loop? –  Mintx Aug 16 '13 at 15:17
    
Weirdly enough, the 85.0511287 value is the max latitude value allowed in EPSG:3857. Because it's the same for all Y values, I would guess a parsing issue with the input data. –  mkennedy Aug 16 '13 at 16:58
    
Maybe its a bad idea to name a variable (proj4) same as a RGeo class (Proj4). Can you report the content of @coords? –  Andre Joost Aug 25 '13 at 7:15
add comment

1 Answer

up vote 3 down vote accepted

The solution--- taken from the creator of RGeo himself, Daniel Azuma, in his rubyconf talk, which is hosted on youtube here: http://www.youtube.com/watch?v=QI0e2jkUbkk, his slides are here: http://speakerdeck.com/u/dazuma/p/getting-down-to-earth-geospatial-analysis-with-rails

in db migration

add_column :model_names, :shape, :spatial, :srid => 2264

in model

class ModelName < ActiveRecord::Base
  proj4 = "+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs"
  GEOFACTORY = RGeo::Geographic.projected_factory(:projection_proj4 => proj4, :projection_srid => 3361)
  set_rgeo_factory_for_column(:shape, GEOFACTORY.projection_factory)
end

to parse:

geo_shape = ModelName::GEOFACTORY.parse_wkt(wkt_string) 
proj_shape = geo_shape.projection

There. proj_shape is the product of geographic coordinates translated to the proper projected coordinates.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.