Zcta sample trouble (Geo-Rails part 8)

[Jose Ernesto Echeverria]

I've been following the RGeo series after listening to Daniel's great presentation at RailsConf. It's really made a difference in the way I'm learning to include GIS elements within my projects.

However, I've come to find something that I'd like to share with you related to Part 8 in the series.

When reading the shapefile as shown in the 1st example (no optimization), I end up with a count of 662 zctas. 

However when reading the shapefile with optimization, I only ended up with 248 zctas

Result is obviously strange as the optimization should have produced a larger number of zctas.

Researching and "prying" into the code I discovered that code as it is for some unexepected reason is not actually fragmenting the zctas but only filtering those that are already below the MAX_SIZE. I discovered that the polygon.intersection(box) isn't producing new polygons for handling the geometry, even though the code looks ok to me.

The following result will illustrate what is happening in my environment, note that it contains a few output instructions for debugging purposes:

1:[zcta=98822,sides=4030][regular4]

[quadrant - min_x:-13449277.23, min_y:6041355.52, max_x:-13414100.39, max_y:6088043.36]

From: /Users/ernestoe/rails/geo_rails_test/db/seeds.rb @ line 90 Object#handle_quadrant:

     90: def handle_quadrant(depth, polygon, min_x, min_y, max_x, max_y, zcta)

     91:   # We do this by creating a rectangle for the box, and computing

     92:   # the intersection with the input polygon. The result could be a

     93:   # polygon, a MultiPolygon, or an empty geometry.

     94:   box = Zcta::FACTORY.polygon(Zcta::FACTORY.linear_ring([

     95:     Zcta::FACTORY.point(min_x, min_y),

     96:     Zcta::FACTORY.point(min_x, max_y),

     97:     Zcta::FACTORY.point(max_x, max_y),

     98:     Zcta::FACTORY.point(max_x, min_y)]))

     99:   #puts box.nil?

    100:   print "\n[quadrant - min_x:#{sprintf("%.2f",min_x)}, min_y:#{sprintf("%.2f",min_y)}, max_x:#{sprintf("%.2f",max_x)}, max_y:#{sprintf("%.2f",max_y)}]"

 => 101:   binding.pry

    102:   handle_geometry(depth, polygon.intersection(box), zcta)

    103: end

[1] pry(main)> box

=> #<RGeo::Geographic::ProjectedPolygonImpl:0x3fcdec8ffd14 "POLYGON ((-13449277.234375233 85.0511287, -13449277.234375233 85.0511287, -13449380.386604048 85.0511287, -13449380.386604048 85.0511287, -13449277.234375233 85.0511287))">

polygon doesn't intersect box, as the newly created box for some unknown reason is producing the same y-coordinate (85.0511287) outside the original range. Notice the y-points of the polygon envelope fall within 6041355 and 6088043.

My guess is something's wrong about projections in my setup, but haven't been able to locate it. I produced the optimization code within the seeds.rb script, don't know whether that is affecting my results. 

If someone that ran the examples can confirm they work, that there is no typo affecting and just in case point me in the right direction, it would be greatly appreciated.

Regards,

--
José Ernesto Echeverría

[ErnestoE]

Kept digging on the same issue and I think I've found the problem.

I changed the code for creating the box:

  box = Zcta::FACTORY.projection_factory.polygon(Zcta::FACTORY.projection_factory.linear_ring([

    Zcta::FACTORY.projection_factory.point(min_x, min_y),

    Zcta::FACTORY.projection_factory.point(min_x, max_y),

    Zcta::FACTORY.projection_factory.point(max_x, max_y),

    Zcta::FACTORY.projection_factory.point(max_x, min_y)]))

Given that we are already using projected coordinates, so we don't need a geographical factory for creating the box but the projection factory.

Now I have 5022 zctas in my database, which seems better than before. 

Can someone please confirm the approach, and if the case, code in the series needs to be corrected.

[Daniel Azuma]

Hi,

I just got around to checking, and it looks like you are correct. The box should be in the projected coordinate system. Whoops!

It's now corrected in the article. Thanks much for investigating and finding this issue.

Daniel