Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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:

Overlaid maps from two servers

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_LIB environment variable, but with no apparent change.
  • The +towgs84=0,0,0,0,0,0,0 part 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.

share|improve this question
    
are you using your own WMS between you client application and the source WMS? –  Shaunak Dec 28 '11 at 12:51
    
Do you mean a proxy WMS server? No - this directly contacts each server. –  David M Dec 28 '11 at 13:23
1  
@DavidM, I checked the conversion in the Esri projection engine and it agrees to 8 decimal places. Beyond that, I don't know. Sorry. –  mkennedy Dec 28 '11 at 20:52
    
@DavidM, Weird speculation: could one be in UTM 55N instead? I checked the test point, and the Y value is about 110km, X is around 300km different. –  mkennedy Dec 28 '11 at 20:57
    
@mkennedy: Thanks for checking! I tried changing the projection to WGS84 UTM 55N (and S) for both maps and for just the problematic one, but they both remain misaligned. (When it's just the problematic one, the coordinates are way off.) When you say you checked the test point and found X and Y off, how did you check it? Ie, even if the ESRI one gives the same as my numbers, what did you do to then discover the reprojected point is wrong? –  David M Dec 29 '11 at 16:52
show 3 more comments

1 Answer 1

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.

share|improve this answer
    
Thanks hmm. I have moved to warping the whole image since I asked this question - not per pixel, but breaking the image into small tiles a few pixels wide. This works great for some projections (UTM maps, for example) but doesn't affect this specific map or projection. –  David M Apr 26 '12 at 12:44
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.