Adding new conversions

[coreation]

Hi!

Great work on this library, it seems to be the only decent one that does coordinate conversion. Can you give me some pointers on what I'd have to do to add more conversions? I'm looking into using UK data, with an EPSG code of 27700, which is not yet supported here. I have no clue of how the math would work to make it go to WGS84 for example. Is there a reference on the web that clarifies the math involved or do I need to seek that out myself reading the specs of the 27700 coordinate system?

[winne27]

Go to http://spatialreference.org

You will find most of epsg codes there. Search your code and than select proj4-code.

Maybe it is agood idea to put this link in README.md.

[judgej]

OSGB uses the OSGB36 datum, which includes the airy ellipsoid. Converting from UK easting/northing values involves a projection transformation first to lat/long, then a datum shift to WGS84.

OSGB uses the tmerc projection, with Proj.4 parameters here:

http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/proj4/

Sorry that's nothing like step-by-step instructions, but should hopefully give you a few clues if you look at some of the example transformations in the test scripts.

[coreation]

Thanks! So transforming the coordinates from an airy ellipsoid to lat/lon, is that a standard operation? And would I then apply a shift of those coordinates using the proj.4 parameters? Or how should I interpret the proj 4 parameters?

[judgej]

Well, sort of ;-)

The UK easting/northing values represent points on a flat map. They map into the Airy datum (and ellipsoid and a centre-of-earth point) using the tmerc projection. That projection will convert between easting/northing, and lat/long coordinates (in both directions). It needs to be given the correct parameters, as defined by the Ordnance Survey, and linked to above.

So that deals with the flat map to/from lat/long conversion.

Once converted to lat/long, that then needs to be datum shifted to the WGS84 datum. Without doing that, your UK locations will be out, sometimes by substantial amounts. A lat/long pair is only meaningful if you know what ellipsoid and earth centre it is based on. The OSGB Airy datum and the WGS84 datum are different, and so the datum needs to be shifted - OSGB lat/long is shifted to or from WHS84 lat/long. The code for doing all that is in the Datrum class.

To shift a lat/long datum, the lat/long is first converted to x/y/z, then those coordinates are rotated and shifted to the WGS84 x/y/z, then those values are converted back to lat/long, and you're done.

That's a bit of background, and I'm explaining it because it is not always obvious what is going on or why you need to do it. What I'm less familiar with is exactly how you do those steps through this package. Perhaps if we can get a worked example put together for your use-case, it can go into the documentation as an example?

[judgej]

I'm going to guess you start like this:

$proj4 = new Proj4php();

$proj_osgb36 = new Proj('+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36', $proj4);

$proj_wgs84 = new Proj('EPSG:4326', $proj4);

$easting = 652709.401;
$northing = 6859290.946;
$point_source = new Point($easting, $northing, $proj_osgb36);

$point_dest = $proj4->transform($proj_wgs84, $point_source);

I've not run that yet, but see how it goes. I'm sure it will be wrong, but we can get something working together. Most of the magic is done in the final transform().

[coreation]

Great stuff @judgej I believe adding a "How to contribute" to the README might do wonders ;) I'll find some time next week to get going with this and see what I bump into. Keep you posted!

[nickolanack]

I thought I'd add a quick note about projection conversions since I just noticed this (I was talking to coreation about a different issue elsewhere).

This is in Proj.php line 202

/**
     * Function: loadFromService
     * Creates the REST URL for loading the definition from a web service and
     * loads it.
     */
    public function loadFromService()
    {
        // Load from web service
        $url = Proj4php::$defsLookupService . '/' . $this->srsAuth . '/' . $this->srsProjNumber . '/proj4/';
        try {
            $this->proj4php->addDef(
                strtoupper($this->srsAuth) . ':' . $this->srsProjNumber,
                $this->proj4php->loadScript($url)
            );
        } catch (Exception $e) {
            $this->defsFailed();
        }
    }

so even if the projection is not 'supported' it should be able to download it, right?.

[julien2512]

Here it is :

$proj4->addDef("EPSG:27700",'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs');

$projOSGB36 = new Proj('EPSG:27700',$proj4);
$projWGS84  = new Proj('EPSG:4326', $proj4);
$pointSrc = new Point(671196.3657,1230275.0454,$projOSGB36);
echo "Source : ".$pointSrc->toShortString()." in OSGB36<br>";
$pointDest = $proj4->transform($projWGS84, $pointSrc);
echo "Conversion : ".$pointDest->toShortString()." in WGS84<br><br>";

[judgej]

@julien2512 Out of curiosity, can you pass the def string direct into new Proj(), or does it always need to be added as a def, with a name, first?

@nickolanack In theory, yes. The def service in this case is http://spatialreference.org/ref/epsg/27700/proj4/

[nickolanack]

@julien2512 btw: I have gotten it to work before without requiring ->addDef(...). although I forked a version of proj4php from a completely different repo than this one...

I mean new Proj('EPSG:27700',$proj4); without first calling ->addDef(...)

[nickolanack]

this is another thing that I tried to do in my own project with some success:

$shpfileProjection=new Proj('PROJCS["NAD_1983_UTM_Zone_17N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-81.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]', $proj4);

[julien2512]

Great ! Let's have a try with OGC WKT defs.

[julien2512]

inline-projections branch on stage ! Travis is working on.

[nickolanack]

ok I'm not sure how this library will work with other wkt definitions, but I had to update the method: parseWKT (in Proj.php 370) to handle parsing definitions recursively.

[nickolanack]

at least to successfully parse the wkt definition I wrote above...

[julien2512]

So it works for at least two definitions... Be sure we will had other tests if we need to !

[coreation]

Aha, great stuff! I'll try the solution @julien2512 gave and let you guys know if it worked. I think this issue can be closed as @nickolanack confirmed it to be working :)

[julien2512]

@nickolanack I miss writing down yours reference in the commit. May I add you on the composer file ? Otherwise you can pull request your Proj.php if you wanted to.

[nickolanack]

You can do that if you want, I don't think I'll PR my other fork, because I made other changes that I don't think you necessarily want, but I will contribute to this project if I can. Thanks

[julien2512]

You're welcome.

I've found the git way to add you on my previous commit. f2aff84ad3c5b0981078477baf087ee290f69430 is now from you.

[coreation]

Hi guys,

I pulled the master into my project so I decided to try custom projections out, starting with a working example, the WGS84 projection.

I'm getting this error: Undefined property: proj4php\Proj::$srsCode on line 447 of proj4php/proj4php/src/Proj.php

The code that triggers this is: new Proj($projCode, $proj4); where $projCode is: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]](OGC WKT code)

As I read the comments of @nickolanack and the merge that happened from his fork I assumed this would be possible, but apparently I'm either doing something wrong, or there's still an issue with initializing the Proj class.

[nickolanack]

oh, I see what is wrong... I think

[coreation]

Mind sharing it? ;)

[nickolanack]

yeah I tried to PR but it failed. can you try to rename the parameter $srsCodeInput to $srsCode and see if it fixes your issue on lines 98 and 110 of Proj.php

[coreation]

That did the trick! Or at least got rid of the error, the transformation itself should also work, I'll leave a comment on this thread if it doesn't. @julien2512 will need to fix it or @nickolanack will have to do a small PR :).

[nickolanack]

I've PR'd my change, which was just that. however the travis ci sever failed testing it so I'll leave it for now, as I'm not sure why it failed

[julien2512]

done, including another test unit for that.

[coreation]

To add more testing to the inline projections, here are a few that fail but of which the EPSG code is already defined. I'm doing this because I want a single way to initiate my projections. So not using the 'EPSG:XXXXX' and 'GEOGCS/PROJCS[....]' at the same time.

These fail:

PROJCS["Belge 1972 / Belgian Lambert 72",GEOGCS["Belge 1972",DATUM["Reseau_National_Belge_1972",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1],AUTHORITY["EPSG","6313"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4313"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",51.16666723333333],PARAMETER["standard_parallel_2",49.8333339],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",4.367486666666666],PARAMETER["false_easting",150000.013],PARAMETER["false_northing",5400088.438],AUTHORITY["EPSG","31370"],AXIS["X",EAST],AXIS["Y",NORTH]]

Error message: Undefined index: Lambert_Conformal_Conic_2SP

GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]]

Error message: geocent:lat out of range:2754.6880650078

This is from the latest pull from master. Edit: haven't tried with others just yet but I think at this point, I'm pretty sure we have some fixing to do. If there's anything I can help with a few pointers would be nice and I'll get started :).

[julien2512]

Good ! Let's start again the inline-projections branch.

First, we need to add your use cases to the php unit (test directory).

For the first one, the error message mean a projection "LCC 2SP" is missing. LCC already exists, it may be similar ?

[coreation]

It is yes, a quick Google search indicated that it had something to do with an offset of 2 parallels. Perhaps that projection calculation is already incorporated in the OGC WKT string?

[nickolanack]

This might help with the lcc issue, although I have no idea if it is the solution, or is even correct... I think that the following parameter definitions are missing for the wkt parsing in Proj.php around line 500

 case 'standard_parallel_1':
   $this->lat1 = $value * Common::D2R;
   break;
 case 'standard_parallel_2':
   $this->lat2 = $value * Common::D2R;
   break;

Otherwise t looks like LCC 2SP as well as LCC 1SP are both handled by lcc.php http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html (defines 2sp and 1sp)

[judgej]

Is there any reason why PHP's deg2rad() and rad2deg() functions can't be used here? I think Common::D2R has just been inherited from JavaScript.

[nickolanack]

they are the same, Common::D2R==deg2rad(1), I don't have a preference other than consistency.

[nickolanack]

@coreation I don't have an error loading the projection:

GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]]

are you getting the error when you transform? if so, what are you transforming to?

[nickolanack]

here is a test that passes:

$proj4           = new Proj4php();

        $projWGS84       = new Proj('EPSG:4326', $proj4);
        $pointSrc = new Point(2.9964931538756,60.863435314163,$projWGS84);
        $newProj=new Proj('GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]]',$proj4);
        $pointNewProj=$proj4->transform($newProj, $pointSrc);
        $pointWGS84=$proj4->transform($projWGS84, $pointNewProj);

        $this->assertEquals(2.9964931538756, $pointWGS84->x, '', 20);
        $this->assertEquals(60.863435314163, $pointWGS84->y, '', 20);

this at least shows that it loads and transforms from, and back to WGS84 (although it does not check the values in the middle...)

[coreation]

@nickolanack I use the triple argument function:

$this->proj4->transform($this->projSrc, $this->projDest, $pointSrc);

Edit: TL;DR a UK geo service contains a wrongly put SHP file/projection code which puts the points next to the equator instead of in the UK. A different dataset from their website however does contain the correct information and is nicely projected using the projection constructor. If the PR of @nickolanack works as well then we're really getting close to closing this issue.

On a smaller note @nickolanack, @julien2512 could it be that this is not working, the srsCode property could not be found, probably because it's protected in the Proj class?

[nickolanack]

yeah your right that is a bug - $srsCode should be made public

[julien2512]

$srsCode is public on inline-projections branch. Any thing else before the master pull request ?

[nickolanack]

I think you should remove this line, I think it is unnecessary. line 110 Proj.php $this->srsCodeInput = $srsCode; it will probably print a warning

[julien2512]

Ok I've finally comment that 110 line. Everything still working. It will be included in the future master commit.

[julien2512]

Well done every one !

The pull request is done.