I don't have a definitive answer, but a step towards a custom CRS using ArcGIS. ArcGIS has a "Local" projection. It actually an ellipsoid-based orthographic projection. The "trans" values were throwing me off until I realized they were false easting/northing values for the local system. This technique isn't useful unless you have control points in both systems, which you provided for nhopton. I would not have gotten this far without him asking for sample points.
Anyway, using the UTM coordinates, unprojected to lat/lon for the center point, and the other parameters, I made a Local CRS. The data still rotated because the rotation value is likely based on the UTM zone, not from geodetic North, so I adjusted it. The points still do not overlay that well (0.2 to 0.6 m) but you can now keep adjusting the parameters to see if you can get a better fit. I just don't have the time right now. Here's the WKT:
PROJCS["canada_local",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["Local"],PARAMETER["False_Easting",2068.553],PARAMETER["False_Northing",4973.397],PARAMETER["Scale_Factor",1.0003933608],PARAMETER["Azimuth",-9.98000000000015],PARAMETER["Longitude_Of_Center",-128.6996571570882],PARAMETER["Latitude_Of_Center",54.0121078922195],UNIT["Meter",1.0]]
Copy the string (as a single line) to a text file and add it to your "Favorites" location. On XP with ArcGIS 10.1, it's
C:\Documents and Settings\login\Application Data\ESRI\Desktop10.1\ArcMap\Coordinate Systems
Other OS or versions will be different.
Now add your UTM or lat/lon data to ArcMap, plus the data in the local system. Do not assign this CRS to the local data. Assign to ArcMap in the data frame properties. You'll see that the reference and local data, almost line up. Now modify the data frame's coordinate system and keep adjusting the parameter values to try to get a better fit. You can use the Apply button on the data frame properties dialog to check how the fit changes.
std disclaimer: I work for Esri.
To visualize the equirectangular projection, I chose this projection string :
+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
and reprojected Natural Earths world boundaries and grid into it using QGIS:
![enter image description here](https://i.stack.imgur.com/BXRAm.png)
The extent of the layers is +/- 20037400 horizontally and half of it vertically. So one degree is about 111km in both directions, what you might call a scale factor to transform directly from angular degrees to projected metres.
The WKT definition as stored in the .prj file is:
PROJCS["Equidistant_Cylindrical",
GEOGCS["GCS_WGS_1984",
DATUM["D_unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Equidistant_Cylindrical"],
PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
If you don't like the results you get with that projection, try a sphere instead of the elliposid. For example EPSG:3786:
+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs
Best Answer
I assume your data are stored inside a PostGIS database. You can use PostGIS functions to do the transformations directly. No need for external libraries.