[gdal-dev] Need some help with projections

[Paul Meems]

Hi list,

I have some data that comes from different sources and are in different
projections.
To use that data I need to reproject them to the same projection.

I'm going to use MapWindow GIS to show the layers.

I've got a XML (not GML) with a lot of info and one of that info is several
coordinates that make up a polygon in epsg:28992 (Dutch RD).
I have a huge data set of shapefiles covering the entire area of The
Netherlands. Those shapefiles are in a custom projection.
And I want to use OpenStreetMap tiles as a background layers.
I have also a GPS receiver connected: WGS84

I'm using this proj4 string for Dutch RD:
+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079
+x_0=155000 +y_0=463000 +ellps=bessel
+towgs84=565.04,49.91,465.84,-0.40939,0.35971,-1.86849,4.0772 +units=m
+no_defs

I'm using this proj4 string for the OSM tiles:
+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m
+no_defs

My WGS84 proj4 string is:
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs

The custom projection used with the shapefiles is to line it up with Google
Earth. That might be working but with OSM the shapefiles are about 10m off.
This is the algorithm used to convert the Dutch RD coordinates is:

DPoint rd2LatLongGE(DPoint rd) //shp_order
{
	// corrections for the rd coordinates for a better match with GE
	const double Xc = 17.78369991;
	const double Yc = 40.18290427;
	
	const double Xx = -0.0000749081;

	const double Xy = -0.0000128123;
	const double Yx = -0.0000064425;
	const double Yy = -0.0000846220;

	return rd2LatLongWgs84(
		DPoint(
			rd.first  + Xc + Xx * rd.first + Xy * rd.second
		,	rd.second + Yc + Yx * rd.first + Yy * rd.second

		)
	);
}


I think the best approach is to convert the shapefiles back to Dutch RD
(their original projection) using a reversed version of the above algorithm.
Can I do that using ogr2ogr? What would the proj4 style project string be?
Don't suggest to get the original shapefiles in DutchRD instead. I can't.
People had to go to court to get this version (
https://rejo.zenger.nl/inzicht/geografische-bestand-3d-shapes all in Dutch).

When I've got the shapefiles in Dutch RD I can check if it is aligning
correctly with my Dutch RD polygon. If it is correct I will convert all data
to epsg:900913 (Google Mercator).

So my main question is what is the proj4 style projection string of the
algorithm and can I use it to reproject the shapefiles.

Thanks,

Paul
The Netherlands
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100714/6eb3f5cc/attachment-0001.html