How do I convert EPSG:26986 coordinate pairs to EPSG:4326 from the MassGIS web service, the Linux command-line, or PHP?

Question

I'm a GIS newbie, normally a PHP web applications developer.

I am trying to map all boundaries of the MA state representative districts. So I just need to convert polygon coordinates from the MassGIS web service's output into Google-friendly lon/lat coordinates (the WGS84 SRS). In GIS parlance I need to execute a "re-projection."

Below is the link to the XML (6MB, careful). The default SRS is EPSG:26986. You can add the GET parameter &srsName=26986 to the end of it and then toy with that number but I've tried plugging in the WGS84 equivalents in EPSG, which folks are saying are 4326 or 3875, and apparently MassGIS does not have those output formats available.

Here is the response XML: http://giswebservices.massgis.state.ma.us/geoserver/wfs?request=GetFeature&version=1.0.0&typeName=massgis:GISDATA.HOUSE2012_POLY

Then I started playing with the Linux tool gdaltransform which seems promising, but when I plugged in over 100 of the original (EPSG:26986) pairs, it only spits back about four pairs (in the correct EPSG:4326 format). So something is wrong there. If tried dumping the coordinate pairs in a text file of source data, but then the gdaltransform complains that it is not a supported file format.

Obviously the easiest thing would be if someone has some experience with the MassGIS web service and can tell me the right parameter to get these coordinates in the format I need them. But another fine solution would be to help me find a simple, quick conversion script or utility (PHP or Linux) that I could integrate into my web app behind the scenes.

Answer 1

Do you know about cs2cs ?

> # usage: cs2cs [proj.source] +to [proj.target]
> cs2cs -r +init=epsg:4326 +to +init=epsg:3857 << EOF
> 45 12
> 45 13
> EOF
1335833.89  5621521.49 0.00
1447153.38  5621521.49 0.00

I translated two coordinates (45ºN,12ºE) and (45ºN,13ºE) to EPSG:3857 (web mercator).

I specified -r (reverse) so to input latitude first (as is defined in EPSG:4326).

The third number in the output triplets is the height Z.

See man cs2cs (APT package is proj-bin here in APT world).

Answer 2

If you put srsName=epsg:4326 instead of just the number you should get a reprojected set of coordinates. This is not standard but relies on a GeoServer specific parameter.

For a more general solution you can use OGR (specifically ogr2ogr) to import the saved file and reproject it. GDAL is designed for rasters rather than vectors which is why you only get the corners.

ogr2ogr -t_srs EPSG:4326 output.shp input.shp