I'm implementing support for WMS maps in my company's software. Servers return a list of projections their maps are available in, but our software only displays data in WGS84 (lat/lon or UTM), and not all servers return the "CRS:84" or "EPSG:4326" CRSes. As a result I need to transform / reproject the data. Some servers return quite unusual ones, and I'm having trouble with one of these: EPSG:3112 (GDA94 / Geoscience Australia Lambert) used by Geoscience Australia's WMS server.
In order to display a map onscreen, I need to convert the WGS84 lat/lon coordinates of the required screen or tile bounding box to the server's CRS, request the map image, and then to draw convert the map's bounding box (in the server's CRS) back to WGS-84 lat/lon. This does not seem to be working quite right:
This image is a portion of the Australian island of Tasmania. The grey coastline is from Demis, in WGS:84 projection. The blue outline and grey graticule is from Geoscience Australia, in EPSG:3112, but (theoretically) reprojected through Proj4 to match the grey coastline. As you can see it's quite inaccurate, to the tune of 50-100 kilometres. If this display were drawing WGS-84 graticules, they would be perfectly vertical and horizontal. Note also that the island is stretched, i.e. it's more than a rotation. This looks to me as though the reprojection is not working at all.
The two Proj4 initialisation strings I'm using are:
- WGS-84 (WMS's CRS:84):
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs - EPSG:3112:
+proj=lcc +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=134 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Round-tripping an individual point from WGS-84 through EPSG:3112 and back yields the original point:
- WGS:84: (145.0, -42.0)
- to EPSG:3112: (931684.430679232, -4764432.89032608)
- and back to WGS-84: (145, -41.9999999999999)
All conversions are done using pj_transform.
I don't know enough about map projections or Proj4 to know where to go next, although I have been learning a lot over the past few days, and I'd really appreciate advice from someone who knows a lot more about this library.
Some ideas I've had are:
- The blue island appears to be still in Lambert projection. Is this likely? It would explain the visual distortion.
- We use a version of the Proj4 DLL that was built years ago, and while it appears you can build various data about datums etc into the DLL, I don't know whether this was done. Proj4 does not flag an error doing the coordinate conversion, though, so at the moment I'm working on the assumption that if the library is happy, it's an external error (i.e. on my part.) I also don't know the reasoning for having our own, years-old build of the DLL anyway, instead of using a prebuilt binary from the Proj4 website. However, even the Proj4 binaries are one version old and don't seem to be built to include everything as binary data.
- Some versions of Proj4 seem to load datum and other information externally. I have tried replacing our DLL and configuring it to use some of these files using the
PROJ_LIBenvironment variable, but with no apparent change. - The
+towgs84=0,0,0,0,0,0,0part of the EPSG:3112 string seems odd. This is from the 'epsg' file in the latest Proj4 download, but could it be wrong?
Perhaps a better formulation of this question is:
- What do I need to have or do to convert between these two CRSes?
Thanks for any help!
I originally posted this question on Stack Overflow, but I've
re-posted it here. I think it's more a GIS question than a
programming question.
Best Answer
I'm not sure if I missed something about your explanation, but normally is not enough to transform only the corners of the bounding box. That only makes match the corners of the image - and not always -, but does not accounts for inner points. You should reproject every pixel in the image to get a perfect match. ¿Maybe using gdalwarp? http://www.gdal.org/gdalwarp.html Another way would be to use a WMS server, for instance Geoserver, as a proxy WMS server to add support for CRS:84 or EPSG:4326 to the services that require it, which is usually known as cascading WMS: http://docs.geoserver.org/latest/en/user/data/wms.html Hope this helps.