Problem with local projections and worldwide data extent #549

Closed
opened this Issue Oct 11, 2011 · 8 comments

None yet

1 participant

Mapnik member

[springmeyer]

``````--- feature_style_processor.hpp 2010-03-10 21:40:34.000000000 +0100
+++ feature_style_processor.hpp.fixed   2010-04-01 23:12:11.000000000 +0200
@@ -126,30 +126,75 @@
proj_transform prj_trans(proj0,proj1);

Envelope<double> layer_ext = lay.envelope();
+
+               // have to determine whether there is an intersection between the map
+               // and the layer extent. since both may have different projections,
+               // one must be reprojected to match the other. this is not something
+               // that always works, e.g. if layer extent = full earth and map uses
+               // a projection only defined for a small area, then projecting the layer
+               // to map's projection will not work. it can also happen the other way
+               // round. thus we try both.
+
+               // first try to back project layer's extent into main map projection
double lx0 = layer_ext.minx();
double ly0 = layer_ext.miny();
double lz0 = 0.0;
double lx1 = layer_ext.maxx();
double ly1 = layer_ext.maxy();
double lz1 = 0.0;
-               // back project layers extent into main map projection
-               prj_trans.backward(lx0,ly0,lz0);
-               prj_trans.backward(lx1,ly1,lz1);
-
-               // if no intersection then nothing to do for layer
-               if ( lx0 > ext.maxx() || lx1 < ext.minx() || ly0 > ext.maxy() || ly1 < ext.miny() )
+               if (prj_trans.backward(lx0,ly0,lz0) && prj_trans.backward(lx1,ly1,lz1))
{
-                  return;
-               }
+                  // if no intersection then nothing to do for layer
+                  if ( lx0 > ext.maxx() || lx1 < ext.minx() || ly0 > ext.maxy() || ly1 < ext.miny() )
+                  {
+#ifdef MAPNIK_DEBUG
+                     std::clog << "no intersection with layer " << lay.name() << " (backward projected)" << std::endl;
+#endif
+                     return;
+                  }

-               // clip query bbox
-               lx0 = std::max(ext.minx(),lx0);
-               ly0 = std::max(ext.miny(),ly0);
-               lx1 = std::min(ext.maxx(),lx1);
-               ly1 = std::min(ext.maxy(),ly1);
+                  // clip query bbox
+                  lx0 = std::max(ext.minx(),lx0);
+                  ly0 = std::max(ext.miny(),ly0);
+                  lx1 = std::min(ext.maxx(),lx1);
+                  ly1 = std::min(ext.maxy(),ly1);

-               prj_trans.forward(lx0,ly0,lz0);
-               prj_trans.forward(lx1,ly1,lz1);
+                  prj_trans.forward(lx0,ly0,lz0);
+                  prj_trans.forward(lx1,ly1,lz1);
+               }
+               // else, if backward projecting layer's extent to map's projection didn't work,
+               // forward project map's extent to layer's projection
+               else
+               {
+                  lx0 = ext.minx();
+                  ly0 = ext.miny();
+                  lz0 = 0.0;
+                  lx1 = ext.maxx();
+                  ly1 = ext.maxy();
+                  lz1 = 0.0;
+                  if (!(prj_trans.forward(lx0, ly0, lz0) && prj_trans.forward(lx1, ly1, lz1)))
+                  {
+                     // this is bad; we can neither project the layer to map nor map to layer.
+#ifdef MAPNIK_DEBUG
+                     std::clog << "cannot render layer " << lay.name() << ": neither forward nor backward projection works" << std::endl;
+#endif
+                     return;
+                  }
+                  if ( lx0 > layer_ext.maxx() || lx1 < layer_ext.minx() || ly0 > layer_ext.maxy() || ly1 < layer_ext.miny() )
+                  {
+#ifdef MAPNIK_DEBUG
+                     std::clog << "no intersection with layer " << lay.name() << " (forward projected)" << std::endl;
+#endif
+                     return;
+                  }
+
+                  // clip query bbox
+                  lx0 = std::max(layer_ext.minx(),lx0);
+                  ly0 = std::max(layer_ext.miny(),ly0);
+                  lx1 = std::min(layer_ext.maxx(),lx1);
+                  ly1 = std::min(layer_ext.maxy(),ly1);
+               }
+
Envelope<double> bbox(lx0,ly0,lx1,ly1);

double resolution = m_.getWidth()/m_.getCurrentExtent().width();
``````
Mapnik member

[olt] The patch is not sufficient. prj_trans.backward can succeed but the return can still be non-sensible. It should try a forward transformation and intersection test if the backward transformation does not intersects.

The result bbox is also to small, some parts are missing. It is not enough to transform the two min and max points, at least the four corners should be used, better ~16 or more points to get a good approximation of the transformed bbox shape.

Mapnik member

[olt] The attached patch adds methods to transform envelopes over the geometry of the envelope and not only the min/max points.
It also tries forward and backward projection like the first patch.

Tested it with global datasets (EPSG:900913 and EPSG:4326) and requests in UTM32N (EPSG:28532) and EPSG:900913.

Mapnik member

[olt] Fixed a bug in the envelope transformation.

Mapnik member

[springmeyer] olt, thanks for the patch. Its needs small improvements (passing Envelope as const&, fixed array length, etc things that will make sure it is fast), but overall the approach looks good. I'll try to find time to test and fix in the next couple of weeks. My goal will be to generate a test suite to make sure that this fixes most/all the other problems listed at BoundsClipping. If you are interested in providing test cases they ideally be a small chunk of data in shapefile or sqlite format + an xml stylesheet that can be added to the nose tests.

Mapnik member

[wohnout] Hello after discussion with springmeyer I'm adding some info here. I'm running debian lenny and I need mapnik as wms (or tilecache) in different SRS than 4326 or 900913. I have installed for this purpose mod_mapnik_wms. OSM file was downloaded from download.geofabrik.de and imported to db using osm2pgsql and style from trunk. I have problem with different SRS (4326 and 900913 is working fine).
http://wms.ccss.cz/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&BBOX=11.025648,46.918073,19.739145,51.689377&SRS=EPSG:4326&WIDTH=1671&HEIGHT=915&LAYERS=OpenStreetMap%20WMS&STYLES=&FORMAT=image/png&DPI=96&TRANSPARENT=true

102067 is <102067> +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813975277778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +towgs84=498.17,136.89,510.08,6.007,4.343,3.831,3.38 +no_defs <>

it's not standard projection (at least at epsg) but we have proof that this definition is working.

I have tried this patch (olt) on 0.7.1, 0.7.2-dev and didn't worked (compile error) and this springmeyer patch workes (it compiles) but no change...

Mapnik member

[NinaG] Hi wohnout,

are you able to solve that problem (compile error)? Your <102067> projection seems to be working.

Mapnik member

[springmeyer] As part of a larger effort to review the tickets at http://trac.mapnik.org/wiki/BoundsClipping I have integrated this approach of first trying to forward the map bbox into the layer srs, then falling back to the original approach of layer back projected into the map srs. Applied in r2784, and will therefore close this ticket. Thanks to frederik for the initial patch and for the other contributions by olt and wohnout.

Note: this fallback should be now rarely needed if the user can supply a 'maximum-extent' for the map in the coords of the map srs. See #506 for more detail. But, a more robust, automatic method is provided in olt's revised patch, which I have finally gotten around to testing. I've moved this portion of the patch (which I did not apply) to a new ticket here: #751, because I feel like it warrants more discussion. Notably it can lead to a significant performance hit for multi-threaded servers. This is not the fault of the patch by rather that proj4 locks for each transform (even if you compile against proj4 trunk and the mapnik mutexes are not used proj4 still locks internally). In tests using paleoserver with 2 threads and proj 4.7 (so using mapniks internal locks) the impact of using 16 points for bbox transformation was 10 r/s (so a tile that could be served at ~30 r/s using the normal 2 point/corner projection and a 'maximum-extent' value dropped to < 20 r/s with the 16 point projection).

closed this Oct 11, 2011
pushed a commit that referenced this issue Oct 11, 2011
 springmeyer `if 'maximum-extent' is provided, clip map query extent to it allowing…` `… one simple and fast way of controlling out of bounds coordinates when projecting into/from layer projection - see also solution at #549` `fd54be4`
pushed a commit that referenced this issue Oct 11, 2011
 springmeyer `first forward map ext to layer srs, falling back to backward projecti…` `…ng layer to map (previous behavior). Also leverage new clip function of box2d class - closes #402, #548, and #549 (and indirectly fixes #308 combined with #506) - see also http://trac.mapnik.org/wiki/BoundsClipping` `69eb5c4`