search

locationtech / geotrellis

GeoTrellis is a geographic data processing engine for high performance applications.
http://geotrellis.io
Other
1.34k stars 361 forks source link

Issue with projecting to EPSG:27700 #1101

Closed lossyrob closed 8 years ago

lossyrob commented 9 years ago

From @mfsny on Gitter

I have my data in the British EPSG:27700 and my leaflet is using WSG84 for passing bbox. I tried following first: 

val src = CRS.fromName("EPSG:3857") 
val dst = CRS.fromName("EPSG:27700") 
val oldx = Extent.fromString("-39135.75848201024,6692214.700423751,0,6731350.458905762") 
val newx = oldx.reproject(src,dst)

but all I get is geotrellis.vector.ExtentRangeError: y: -1.3747887674479786E7 to -1.374961923978705E7.

So I checked the proj4 string in /nad/epsg and it looks like the value (<27700> +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs <>) is not correct. I found some other defs online and this one worked the best: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs. But it has still some southwards. I get 514302.4577355564,188783.730894182,537991.4618278733,213686.72817825677 when I reproject using geotrellis. http://cs2cs.mygeodata.eu/ says it should be : 514771.247816,167974.826075,538560.963803,192915.540317 .

dwins commented 8 years ago

Looks like proj4j has no implementation of grid shifting: https://github.com/geotrellis/geotrellis/blob/70a708a/proj4/src/main/java/org/osgeo/proj4j/BasicCoordinateTransform.java#L162 . A debug session with proj.4 indicates this would be required for correct reprojection in this case. I verified that replacing the grid shift implementation (pj_apply_gridshift2) in proj.4 with a no-op produces identical results to what we're getting from proj4j, but it looks like porting that code over to proj4j is going to be a decent-sized undertaking.

dwins commented 8 years ago

After reviewing this on my branch, I believe that merging PR #1306 will resolve the issue. Note that in reprojecting envelopes GeoTrellis reprojects all 4 corners and takes the envelope of the resulting polygon, so only reprojecting the lower left and upper right corners as you have gives a different result. However I see results in line with what you expect when I do a point-wise reprojection:

val src = CRS.fromName("EPSG:3857") 
val dst = CRS.fromName("EPSG:27700") 
val dstAlt = CRS.fromString("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157
val oldx = List[Point](
  (-39135.75848201024,6692214.700423751),
  (0d,6731350.458905762)) 
val newx = oldx.map(_.reproject(src,dst))
val altx = oldx.map(_.reproject(src, dstAlt))
println(newx)
println(altx)

Produces 

List(POINT (514771.25069103093 167974.82616902422), POINT (538560.9667028334 192915.54038821382))
List(POINT (514771.2478158922 167974.82607492874), POINT (538560.9638029323 192915.54031725624))