Mismatch with vector projection

[T. D.]

Hello everyone 

After working a while with tilestache and searching a bit over the net, I couldn't find a solution to my problem. I am using a pgsql/postgis db imported from OSM with osm2pgsql with Sperical Mercator projection. I want to create 2 layers, one from the osm_planet_line with mapnik and png tiles, and one from osm_planet_point with geojson tiles. I have successfully connected the db with tilestache and created the layers, I have set the projection to spherical Mercator, but there seems to be a mismatch between the layers. I centered the map on a specific point (expressed in latitude and longitude) and compared the result that was produced by mapnik with the one from the official OSM and they match. So I believe I have done something wrong with the configuration of the vector layer.

Here is my cfg:

{

  "cache": {

    "name": "Disk",

    "path": "/tmp/tiles/",

    "dirs": "portable",

    "gzip": [],

    "logging": "debug"

  },

  "layers": {

    "osm_point":

    {

      "provider": {"name": "vector", "driver": "postgis",

                   "parameters": {"dbname": "gis", "user": "postgres", "password": "12345",

                                  "table": "planet_osm_point"}},

      "projection": "spherical mercator",

      "allowed origin": "*",

      "bounds": { "low":9 }

    },

    "osm_ori_line":

    {

        "provider": {"name": "mapnik", "mapfile": "/home/ubuntu/mapnik2.xml"}

    }

  }

}

In case you need it, here is the mapnik configuration:

<?xml version="1.0" encoding="utf-8"?>

<!DOCTYPE Map [

<!ENTITY % entities SYSTEM "/home/ubuntu/mapnik/inc/entities.xml.inc">

%entities;

]>

<Map background-color="#b5d0d0" srs="&srs900913;" minimum-version="2.0.0">

  <Style name="MyStyle">

    <Rule>

      <LineSymbolizer stroke="rgb(10%,10%,10%)" stroke-width="0.2" />

    </Rule>

  </Style>

<Layer name="osm_ori_line" status="on" srs="&srs900913;">

    <StyleName>MyStyle</StyleName>

    <Datasource>

        <Parameter name="type">postgis</Parameter>

        <Parameter name="host">localhost</Parameter>

        <Parameter name="dbname">gis</Parameter>

        <Parameter name="user">postgres</Parameter>

        <Parameter name="password">12345</Parameter>

      <Parameter name="table">

      (select way

       from  planet_osm_polygon

        UNION

        select way

        from planet_osm_roads

      ) as MyStyle

      </Parameter>

      &datasource-settings;

    </Datasource>

</Layer>

</Map>

Any ideas would be a great help.

Thank you for your time

Th.D.

[Michal Migurski]

Hi TD,

Can you post an image of the mismatch, and a sample of your dataset?

-mike.

> --
> You received this message because you are subscribed to the Google Groups "Tilestache" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to tilestache+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

----------------------------------------------------------------
michal migurski- contact info and pgp key:
sf/ca http://mike.teczno.com/contact.html

[T. D.]

Sure here is a screenshot http://tinypic.com/r/3093os1/6 . The missmatch is about 50km and i think its only on Longitude

i am attaching a txt file with the results of the query: select * from planet_osm_roads limit 3;

also here is the projection of data

gis=# select getsrid(way) from planet_osm_roads limit 2 ;

 getsrid

---------

  900913

  900913

(2 rows)

Do you need anything else from the dataset?

[Michal Migurski]

What are you using for your Mapnik SRS?

This kind of north-south error is something I've often seen with an incorrectly-defined spherical mercator projection.

The first row in your attached database extract is this road:
http://www.openstreetmap.org/browse/way/95388137

And the extent of the first shape in the extract matches:
http://www.openstreetmap.org/?box=yes&bbox=19.701742%2C39.678700%2C19.728045%2C39.671799

What's in your /home/ubuntu/mapnik/inc/entities.xml.inc ?

-mike.

> <out.txt>

[T. D.]

The srs that i am using is dedined in /home/ubuntu/mapnik/inc/entities.xml.inc:

<!ENTITY srs900913 "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over">
<!ENTITY srsmercator "+proj=merc +datum=WGS84 +over">
<!ENTITY srs4326 "+init=epsg:4326">
<!ENTITY maxscale_zoom0 "<MaxScaleDenominator>250000000000</MaxScaleDenominator>">
..............
<!ENTITY minscale_zoom18 "">

here is the definiton of google projection in /usr/share/proj/epsg in case you need it:

# Google proj

<900913> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over<>

[Michal Migurski]

Well that looks kosher. Any chance the thing you're registering against is wrong?

-mike.

> <entities.xml.inc.txt>

[T. D.]

I am sorry can you explain a bit more what you want me to check?

The thing that troubles me is that both layers use the same data, so it has to be one of the projection rules. How can i check if the vector layer gets the right settings?

[Michal Migurski]

To check if the vector layer has the right *output*, you can compare it against this:
http://teczno.com/squares/

You'll see the correct values for both lat/lon and mercator there, and you can zoom in quite far to match them up to specific points such as intersections or POIs. I'm almost 100% certain you have a datum mismatch, and it usually takes me some trial and error to figure those out.

-mike.

[T. D.]

Hi Mike

Thanks for the fast responses 

I used the site and firebug to test the tiles with their data. I centered to a point and got the tile id that was at that point. I then checked with firebug to see what tiles are requested. Both layers (vector and mapnik) request the same set of tiles, with the same IDs. The mapnik layer centers at the same point, as the link you showed me did. I checked the tile with the same ID from the vector layer and it has different data in it.

I also checked the code from Tilestache to see where the rule for the spherical projection exists. I think I spotted it under Tilestache/Geography.py. I compared the rules between them and it seems that the mapnik rule is missing the +wktext option. I added it, cleared the cache but nothing changed. I also updated from 1.41 to 1.44 still no luck…

Thanks again

T.D.

[T. D.]

Hi all again

I calculated that the missmatch is 0.18667 points in latitude. I know i ask a lot, but does anyone knows a way to "fix" this offset? Should i change the projection srs in a way that it will produce the tiles alligned?

If anyone has an idea on how to fix the datum mismatch please say so. 

Thanks again for the help, Sorry for beeing so persistant on this, it is one thing i need to complete my thesis.

T.D.

[T. D.]

Hello all, please ignore the previous message, its not the solution. I just jumped to conclutions but they were wrong. I deleted it also from the threat.

I couldn't find a way to fix the datum, so i thought to fix it manually throught the data. 

So i cloned the whole database and created a program that connects to the database uses the ST_AsText to take the geometry data, then parse it, convert the Mercator meters to coordinates, add the offset, convert back to Mercator meters, then to geometry with ST_GeomFromText and then store it to another table (named roads2).

Then i created an index on way column but with USING gist option, and inserted to geometry_columns the new geometry column. 

After that I managed to allign the layers. At start i thought that everything was ok, put after panning a lot i figured that a slight missmatch started again, which means that the offset was not a constant value.

I feel that i am pretty close into finding a solution. Any suggestions would be great!

Thanks

T.D.

[Andrew Harvey]

I'm not exactly sure if this will help you, but just letting you know
that I have previously encountered a vertical discrepancy between a
Mapnik and Vector layer in TileStache (with the Vector one being wrong).

I'm not exactly sure what I did to fix it, but I believe it was related
to the projection I was using for the data in PostGIS (so the Vector
layer was pulling from PostGIS).

I believe in the end following this advice from the TileStache
documentation seemed to help.

If you are using PostGIS and spherical mercator a.k.a. SRID 900913,
you can save yourself a world of trouble by using this definition:

http://github.com/straup/postgis-tools/raw/master/spatial_ref_900913-8.3.sql

[T. D.]

Hello Andrew 

and thank you for your help and your respnonce, you are a lifesaver. It actually fixed it. I must have missed that information from the documentation. I had created a 900913 rule following the instructions of postgis installation, but it was a bit different. So i deleted the old rule and created the new one.

Anyway it works perfect now, thanks again everyone for the support, especially Mike and Andrew.

T.D.

[Michal Migurski]

I'm glad this worked for you, and thanks Andrew for picking up the projection issue.

So you had the correct geometry in the database which we verified with that Greek road, and the correct projection in your map which we saw, but the *incorrect* projection in your spatial references table in PostGIS causing Mapnik to reproject things to the wrong place, fixed by the 900913 below.

-mike.

[Michal Migurski]

I think it should be possible to add a test to the vector provider for this. Any chance you could share the 900913 projection you had in your DB from before, or the version of PostGIS/Postgres you used and how you installed it?

-mike.

[Jerry Thomas]

May I just add that while TileStache in concept is great, it seems the entry pains to get it configured are damn near a nightmare for nearly everyone that chooses to utilize it.  What I would like to see is someone with some mad programming skills transform this TileStache thing from a kluge of python scripts and dependencies into a single self-contained executable, with builds for each platform.  Vector based base tools are the future of mapping technologies and I think it is about time that great strides are made towards transforming this into a more mature program with less entry learning curve configuration requirements, and if that meant that it would need to become a paid commercial app, then so be it.  If I had the skills to make this transformation then I would do it.

[Michal Migurski]

Say more, Jerry!

In this case, the problem was found to exist in PostGIS with the mismatched projection, so I'm not sure a packaged Tilestache would have made any difference.

I'm curious, though, what kind of packaging would be most desirable. I imagine that Tilestache as commercial software would miss the mark, but a solid and reliable vector tile service would be a lot closer, something with baked-in and up-to-date OSM data included and ready for use. If vector tiles are the future, do you envision a greater need for one true source of everybody's vector tiles or software for hacking out custom vector tiles?

-mike.

[T. D.]

Yes mike this was the record i had. I installed the 9.2 postegres and  postgis 1.5.3

here is the SELECT PostGIS_full_version(); result:

POSTGIS="1.5.3" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.8" USE_STATS 

as for the record in spatial_ref_sys:

900913 | spatialreferencing.org |    900913 | PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]] | +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs

To be honest i am not 100% sure how i did the installation, cause it was back on September, but if I recall corectlly both postgresql and postgis were installed from packages 

and GDAL and GEOS were installed from source in order to enable the python hooks(there was a guy who listed how to do this in the group).

[T. D.]

If I may share an opinion here, I recently started working with tilestache and i have done great things that would take me months to achieve with any other tools. I started from scratch and probably the most difficult part of getting started is that there isn't any basic guides. The documentation is great, you can find everything in it. But when you start from 0 like i did, you don't know where to start. 

I think that a script that checks requirements, or an installation script that will use specific packages in order to make a basic installation work, would be a great help. 

Btw what i find intresting is that the projection rule in the database was not used by mapnik (mapnik uses its own rules), and generated the tiles right. I think it would be helpfull if there was an option in the layer cfg to use a "fixed" projection rule to avoid future errors like this.

T.D.

[Andrew Harvey]

On 04/03/13 04:30, Michal Migurski wrote:
> I think it should be possible to add a test to the vector provider
> for this. Any chance you could share the 900913 projection you had in
> your DB from before, or the version of PostGIS/Postgres you used and
> how you installed it?

I could be wrong, but I believe I ran into this problem by loading the
geometries into PostGIS using a new PostGIS enabled DB with the
spatial_ref_sys.sql file loaded then using ogr2ogr -t_srs 'EPSG:3857'.

On 04/03/13 04:36, Jerry Thomas wrote:
> it seems
> the entry pains to get it configured are damn near a nightmare for
> nearly everyone that chooses to utilize it. What I would like to see
> is someone with some mad programming skills transform this TileStache
> thing from a kluge of python scripts and dependencies into a single
> self-contained executable

I find it pretty darn simple myself. On Debian unstable you can install
with,

# apt-get install tilestache

Then find the sample configuration at

/usr/share/doc/tilestache/examples/tilestache.cfg

The excellent documentation at

/usr/share/doc/tilestache/html/index.html

Then run the application using either

$ tilestache-server
$ gunicorn "TileStache:WSGITileServer('/path/to/tilestache.cfg')"

or some other method as described in the documentation.

> Vector
> based base tools are the future of mapping technologies

Indeed for myself I wish that this patch would make it into upstream,

https://github.com/andrewharvey/TileStache/commit/1869a509b190a69b22b40d5ba1cbf2b34595bba6

For me it is the killer feature as it allows web map tiled vector data
generalised to your zoom level.

I don't understand why it can't me merged upstream and it would mean I
don't need to maintain my own builds for this feature and it would lead
to further widespread adoption of zoom dependent generalised vector tiles.

[Michal Migurski]

On Mar 4, 2013, at 12:49 AM, Andrew Harvey wrote:

> On 04/03/13 04:30, Michal Migurski wrote:
>> I think it should be possible to add a test to the vector provider
>> for this. Any chance you could share the 900913 projection you had in
>> your DB from before, or the version of PostGIS/Postgres you used and
>> how you installed it?
>
> I could be wrong, but I believe I ran into this problem by loading the
> geometries into PostGIS using a new PostGIS enabled DB with the
> spatial_ref_sys.sql file loaded then using ogr2ogr -t_srs 'EPSG:3857'.

I won't claim this is anything other than plain superstition, but the non-900913 competing 38XX definitions of spherical mercator have consistently led me into trouble and I avoid them. We're coming up on eight years of this incredibly simple, obvious projection being "hard."

>> Vector
>> based base tools are the future of mapping technologies
>
> Indeed for myself I wish that this patch would make it into upstream,
>
> https://github.com/andrewharvey/TileStache/commit/1869a509b190a69b22b40d5ba1cbf2b34595bba6
>
> For me it is the killer feature as it allows web map tiled vector data
> generalised to your zoom level.
>
> I don't understand why it can't me merged upstream and it would mean I
> don't need to maintain my own builds for this feature and it would lead
> to further widespread adoption of zoom dependent generalised vector tiles.

Missed this patch, sorry. I've implemented similar things myself a number of times, I'll prioritize getting it merged.

-mike.

[Paul Norman]

> From: tiles...@googlegroups.com [mailto:tiles...@googlegroups.com]
> On Behalf Of Michal Migurski
> Sent: Monday, March 04, 2013 10:28 AM
> Subject: Re: [tilestache] Mismatch with vector projection

>
> On Mar 4, 2013, at 12:49 AM, Andrew Harvey wrote:
>
> > On 04/03/13 04:30, Michal Migurski wrote:
> >> I think it should be possible to add a test to the vector provider
> >> for this. Any chance you could share the 900913 projection you had in
> >> your DB from before, or the version of PostGIS/Postgres you used and
> >> how you installed it?
> >
> > I could be wrong, but I believe I ran into this problem by loading the
> > geometries into PostGIS using a new PostGIS enabled DB with the
> > spatial_ref_sys.sql file loaded then using ogr2ogr -t_srs 'EPSG:3857'.
>
> I won't claim this is anything other than plain superstition, but the
> non-900913 competing 38XX definitions of spherical mercator have
> consistently led me into trouble and I avoid them. We're coming up on
> eight years of this incredibly simple, obvious projection being "hard."

Having just gone through a set of imagery conversions, the key to getting it
to work it to stick to one projection all the way through the toolchain.
I've had more luck with EPSG:3857 in new software than the old definitions.

[erik...@gmail.com]

Thank you Andrew for the help below! 

This saved me a huge amount of time/hair. 

My vector layer was also way "too much" north than expected. Changing the spatial reference definition in the postgis database (spatial_ref_sys) to the one you supplied solved this issue. 

Andrew/mike: Do you know if it is the "Mercator_2SP" vs "Mercator_1SP" that is the the difference-maker?

mike: Maybe you could add a paragraph on how to check this in your postgis installation in the docs? Maybe even on the "/doc" page so it is more visible?

I did something like: "SELECT srtext FROM spatial_ref_sys WHERE srid = 900913;"

There I checked if I had "Mercator_2SP" or not. If not, delete row ("DELETE FROM spatial_ref_sys WHERE srid = 900913;") and add new row from file below.

Thank you again!

[Max Demars]

 You can read a nice post on that topic here: http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/