Commits

Simon Kokkendorff committed e45111f

Expanded on the proj4 to mlb translation and added tests in testsuite (TestProj.py,AutoTest.py).
The implementation is still shaky regarding datum shifts - sadly this seems to be unavoidable...
TODO: implement the easier way around - from mlb to proj4.

Comments (0)

Files changed (7)

TR_INC/parse_def_file.h

 	char mlb[MLBLNG];
 	int no;
 	char p_datum[MLBLNG];
-	char ellipsoid[MLBLNG];
+	char ellipsoid[MLBLNG]; /*really a reference to a def_grs struct */
 	int imit;
 	int type;
 	int p_no;

TR_INC/proj4_to_mlb.h

  struct proj4_entry{
 	char proj[64];
 	char datum[64];
+	char ellps[64];
 	char units[16];
+	double towgs84[7];
+	int n_params;
 	int zone;
 	double x_0;
 	double y_0;
 	double k;
  };
  
+ 
+	 
+	 
  typedef struct proj4_entry proj4_entry;
+
  
  /*Remember to free the returned pointer*/
  proj4_entry *parse_proj4(char *proj4_text);

TR_SRC/proj4_to_mlb.c

 /*
- * Copyright (c) 2011, National Survey and Cadastre, Denmark
+ * Copyright (c) 2012, National Survey and Cadastre, Denmark
  * (Kort- og Matrikelstyrelsen), kms@kms.dk
  * 
  * Permission to use, copy, modify, and/or distribute this software for any
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
+#include "geo_lab.h"
 #include "proj4_to_mlb.h"
 #include "strong.h"
 #include "trlib_api.h"
 #include "lord.h"
-/* A sketcy implementation of proj4 to mlb translation. 
+#include "parse_def_file.h"
+/* A sketchy implementation of proj4 to mlb translation. 
 * lots of hardcoded stuff, which should perhaps be made more visible....
 */
 
 
+
+struct translation_entry{
+	char *proj4_name;
+	char *kms_name;
+};
+
+struct towgs84_translation{
+	 char *kms_dtm;
+	 char *proj4_ellps;
+	 int min_params;
+	 double towgs84[7];
+ };
+
+ /* not used at the moment...
+static struct proj4_to_kms_ellps 
+ELLIPSOID_TRANSLATIONS[]={
+{"intl","Hayford"},
+{"bessel","Bessel"},
+{NULL,NULL}
+};
+*/
  
+/* the datums that we want to guess from proj4-ellipsoids */
+static struct translation_entry
+ELLIPSOID_TO_DATUM[]={
+{"GRS80","etrs89"}, /*could be gr96 as well */
+{"WGS84","wgs84"},
+{NULL,NULL}
+};
+
+/*translations of proj4 datums to kms datums*/
+static struct translation_entry
+PROJ4_TO_KMS_DATUM[]={
+{"WGS84","wgs84"},
+{"NAD83","nad83"},
+/*perhaps we should add some more ourselves, even if they are not supported by the library proj, but only by the language proj?? */
+{NULL,NULL}
+};
+
+/* A few datums that we want to recognize from ellipsoid and common proj4 datum shift parameters */
+static struct towgs84_translation
+COMMON_PROJ4_DEFINITIONS[]={
+{"ed50","intl",3,{-87,-98,-121,0,0,0,0}},
+{"etrs89","GRS80",3,{0,0,0,0,0,0,0}},   /*and this could be gr96 as well! */
+{NULL,NULL,0,{0,0,0,0,0,0,0}}
+};
+
+	
+
+/*Function that parses proj4 tokens into a struct*/
+
  proj4_entry *parse_proj4(char *proj4_text){
 	 char **tokens=NULL, *no_spaces;
 	 int item_count,i,pair;
 	 /*set default values*/
 	 proj_entry->proj[0]='\0';
 	 proj_entry->datum[0]='\0';
+	 proj_entry->ellps[0]='\0';
 	 strcpy(proj_entry->units,"m");
 	 proj_entry->x_0=0;
 	 proj_entry->y_0=0;
 	 proj_entry->k=1.0;
+	 proj_entry->n_params=0;
 	 proj_entry->lon_0=-9999;
 	 proj_entry->lat_0=-9999;
 	 proj_entry->zone=-9999;
 		else if (!strcmp(key_val[0],"datum"))
 			strncpy(proj_entry->datum,key_val[1],64);
 		else if (!strcmp(key_val[0],"ellps"))
-			strncpy(proj_entry->datum,key_val[1],64);
+			strncpy(proj_entry->ellps,key_val[1],64);
 		else if (!strcmp(key_val[0],"units"))
 			strncpy(proj_entry->units,key_val[1],16);
 		else if (!strcmp(key_val[0],"zone"))
 			proj_entry->lon_0=atof(key_val[1]);
 		else if (!strcmp(key_val[0],"k_0") || !strcmp(key_val[0],"k"))
 			proj_entry->k=atof(key_val[1]);
+		else if (!strcmp(key_val[0],"towgs84")){
+			char **params;
+			int n_params,j;
+			params=get_items(key_val[1],",",&n_params);
+			for(j=0; j<n_params; j++)
+				proj_entry->towgs84[j]=atof(params[j]);
+			free(params);
+			proj_entry->n_params=n_params;
+		}
 		free(key_val);
 	}
 	free(tokens);
 	return proj_entry;
 }
 
+/*Guess a datum from proj4 data 
+* This is a mess if the +datum token is not give explicitely, since there does not seem to universal agreement on the +towgs84 7-parameter transformation for many datums!
+* The proj4 policy seems also to have been not to include the datum-shift in the above case 
+* Thus there is no bullet proof method of setting a KMS datum from a proj4 string, unless the datum is explicitly given with the +datum token.
+*  E.g. there is no way of distinguishing etrs89 and gr96 in proj4 currently. (towgs both zeroes...)
+*  A *few* shaky defaults will be built in here however.... 
+*/
+static int guess_proj4_datum(proj4_entry *proj_entry, char *datum){
+	/*if no datum - look for towgs84 params or guesssssss*/
+	datum[0]='\0';
+	/*if proj4 datum token is not set, we will conduct a search for a KMS-match */
+	if (proj_entry->datum[0]=='\0'){
+		/*if no datum shift parameters are set, guess datum from ellipsoid name from a list of a few defaults - shaky! */
+		if (proj_entry->n_params==0){
+			lord_warning(1,"Proj4: no proper datum (or datum shift) definition. Guessing datum from ellipsoid.");
+			/* if ellipsoid is not set, default to wgs84*/
+			if (proj_entry->ellps[0]=='\0'){
+				lord_warning(1,LORD("Proj4: proj and ellps keys both undefined, setting datum to WGS84"));
+				strcpy(datum,"wgs84");
+			}
+			else{/*look for a default translation of ellipsoid name to datum */
+				struct translation_entry *current_try=ELLIPSOID_TO_DATUM;
+				while (current_try->proj4_name!=NULL){
+					if (!strcmp(current_try->proj4_name,proj_entry->ellps)){
+						strcpy(datum,current_try->kms_name);
+						break;
+					}
+					current_try++;
+				}
+			}
+			
+		}
+		/* else shift parameters are set. Then see if datum shift looks like something we know */
+		else{ 
+			
+			struct towgs84_translation *current_try=COMMON_PROJ4_DEFINITIONS;
+			int i;
+			while (current_try->kms_dtm!=NULL){
+				if (!strcmp(proj_entry->ellps,current_try->proj4_ellps) && proj_entry->n_params>=current_try->min_params){
+					int nok=0;
+					for (i=0; i<proj_entry->n_params && i<7; i++)
+						nok+=(fabs(proj_entry->towgs84[i]-current_try->towgs84[i])<1e-7);
+					if (nok==proj_entry->n_params){
+						strcpy(datum,current_try->kms_dtm);
+						break;
+					} /*end found datum*/
+				} /*end start checking ccommon proj4 def*/
+				current_try++;
+			} /*end loop over common proj4 definitions */
+			/*TODO: perhaps loop through kms-definitions in DEF_DATA to find a match from parameters...*/
+		}
+			
+			
+		
+	} /*end datum token not given */
+	else{ /*datum token is given and we look for a translation */
+	/*TODO: add more of predefined proj4-datums*/
+		struct translation_entry *current_try=PROJ4_TO_KMS_DATUM;
+		while (current_try->proj4_name!=NULL){
+			if (!strcmp(current_try->proj4_name,proj_entry->datum)){
+				strcpy(datum,current_try->kms_name);
+				break;
+			}
+			current_try++;
+		} /*end loop over translations */
+		if (strlen(datum)==0){ /*if datum not supported */
+			lord_error(TR_LABEL_ERROR,"Proj4-datum: %s, not supported.",proj_entry->datum);
+			return TR_LABEL_ERROR;
+		}
+	}/* end else datum is given */
+	return TR_OK;
+}
+/*the actual conversion of proj4 text to mlb */
+
 int proj4_to_mlb(char *proj4_text, char *mlb){
-	int is_transverse_mercator=0;
+	extern def_data *DEF_DATA;
+	int is_transverse_mercator=0,rc;
+	char proj[32],datum[32];
 	proj4_entry *proj_entry=parse_proj4(proj4_text);
+	proj[0]='\0';
+	datum[0]='\0';
 	if (proj_entry==NULL){
 		lord_error(TR_LABEL_ERROR,"Unable to parse proj4 string: %s",proj4_text);
 		return TR_LABEL_ERROR;
 	/*set projection*/
 	
 	if (proj_entry->proj[0]=='\0'){
-		lord_error(TR_LABEL_ERROR,"Proj4: projection not set!");
+		lord_error(TR_LABEL_ERROR,LORD("Proj4: projection not set!"));
 		free(proj_entry);
 		return TR_LABEL_ERROR;
 	}
 	if (!strcmp(proj_entry->proj,"utm")){
-		strcpy(mlb,"utm");
-		if (proj_entry->zone<0)
-			return 1;
-		sprintf(mlb+3,"%2d",proj_entry->zone);
-		mlb+=5;
+		if (proj_entry->zone<0 || proj_entry->zone>60){
+			lord_error(TR_LABEL_ERROR,LORD("Proj4: utm-zone not properly defined"));
+			free(proj_entry);
+			return TR_LABEL_ERROR;
+		}
+		sprintf(proj,"utm%02d",proj_entry->zone);
+		
 	}
 	else if (!strcmp(proj_entry->proj,"latlong") || !(strcmp(proj_entry->proj,"longlat"))){
-		strcpy(mlb,"geo");
+		strcpy(proj,"geo");
 		mlb+=3;
 	}
 	else if (!strcmp(proj_entry->proj,"tmerc") || !(strcmp(proj_entry->proj,"etmerc"))){
-		strcpy(mlb,"itm");
 		is_transverse_mercator=1;
-		mlb+=3;
+		strcpy(proj,"itm");
 	}
 	else{
 		strcpy(mlb,"err");
-		lord_error(TR_LABEL_ERROR,"Proj4-projection: %s, not supported.",proj_entry->proj);
+		lord_error(TR_LABEL_ERROR,LORD("Proj4-projection: %s, not supported."),proj_entry->proj);
 		free(proj_entry);
 		return TR_LABEL_ERROR;
 	}
 	/*set datum*/
-	/*if no datum - set default*/
-	if (proj_entry->datum[0]=='\0'){
-		lord_warning(1,"Proj4-string, no datum, setting datum to WGS84.");
-		strcpy(proj_entry->datum,"WGS84");
-	}
-	if (!strcmp(proj_entry->datum,"WGS84"))
-		mlb+=sprintf(mlb,"_wgs84");
-	else if (!strcmp(proj_entry->datum,"GRS80"))
-		mlb+=sprintf(mlb,"_etrs89");
-	else if (!strcmp(proj_entry->datum,"ED50"))
-		mlb+=sprintf(mlb,"_ed50");
-	else{
-		lord_error(TR_LABEL_ERROR,"Proj4-datum: %s, not supported.",proj_entry->datum);
+	rc=guess_proj4_datum(proj_entry,datum);
+	
+	if (strlen(datum)==0){
+		lord_error(TR_LABEL_ERROR,"Proj4: Unable to recognize datum...");
+		free(proj_entry);
 		return TR_LABEL_ERROR;
 	}
 	/*put extra parameters*/
-	if (is_transverse_mercator)
-		sprintf(mlb,"  %.2fdg  %.2f%s  %.2fdg  %.2f%s  %.8f",proj_entry->lat_0,proj_entry->y_0,proj_entry->units,
-	proj_entry->lon_0,proj_entry->x_0,proj_entry->units,proj_entry->k);
+	if (is_transverse_mercator){
+		/*search for definition*/
+		def_projection *prj=NULL,*this_prj;
+		int i,params_ok;
+		if (DEF_DATA==NULL)
+			lord_warning(1,LORD("definition data not initialised - unable to search for minilabel."));
+		else{
+			for(i=0;i<DEF_DATA->n_prj; i++){
+				params_ok=0;
+				this_prj=DEF_DATA->projections+i;
+				if (this_prj->q_par==5 && (!strcmp(this_prj->native_proj,"itm") || !(strcmp(this_prj->native_proj,"itmi")))) {
+					/*test parameters*/
+					double scale=1;
+					if (!strcmp(proj_entry->units,"km"))
+						scale=1000;
+					params_ok+=(proj_entry->k==this_prj->native_params[4]);
+					params_ok+=(fabs(proj_entry->lat_0-this_prj->native_params[0]*DOUT)<0.0000001);
+					params_ok+=(fabs(proj_entry->y_0*scale-this_prj->native_params[1])<0.000001);
+					params_ok+=(fabs(proj_entry->lon_0-this_prj->native_params[2]*DOUT)<0.0000001);
+					params_ok+=(fabs(proj_entry->x_0*scale-this_prj->native_params[3])<0.000001);
+					if (params_ok==5){
+						prj=this_prj;
+						break;
+					}
+					
+				}
+			}
+		} /*end else*/
+		if (prj!=NULL && !strcmp(datum,prj->p_datum))
+			strcpy(mlb,prj->mlb);
+		else{
+			sprintf(mlb,"itm_%s  %.2fdg  %.2f%s  %.2fdg  %.2f%s  %.8f",datum,proj_entry->lat_0,proj_entry->y_0,proj_entry->units,
+			proj_entry->lon_0,proj_entry->x_0,proj_entry->units,proj_entry->k);
+			if (prj!=NULL)
+				lord_info(1,"Proj4: %s found semi-matching mlb %s, however no datum match.",proj4_text,prj->mlb);
+		}
+	}
+	else{
+		sprintf(mlb,"%s_%s",proj,datum);
+	}
 	free(proj_entry);
 	return TR_OK;
 }

python/AutoTest.py

 import TestFiles
 import RandomTests
 import ThreadTest
+import TestProj
 import time,os,sys,platform
 TEST_DATA=os.path.join(os.path.dirname(__file__),os.path.join("..","testdata"))
 JOB_DEF=os.path.join(TEST_DATA,"job_new.txt")
 		nerr+=1
 		print("Error:\n%s" %repr(msg))
 	print LS
+	try:
+		nerr+=TestProj.main(["TestProj.py"])
+	except Exception,msg:
+		nerr+=1
+		print("Error:\n%s" %repr(msg))
+	print LS
 	#forbid unsafe transformations#
 	TrLib.SetThreadMode(True)
 	try:

python/TestProj.py

+#!/usr/bin/python
+#######################################
+##Skeleton of a test suite for trlib using python 'bindings'
+## This script performs some translations of proj4 definitions to kms-minilabels
+## simlk, oct. 2012
+#######################################
+import sys,os,time
+import TrLib
+TESTS=["+proj=utm +zone=32 +units=m +ellps=GRS80 +nodefs",
+"+proj=utm +zone=32 +units=m +datum=WGS84 +nodefs",
+"+proj=latlong +ellps=GRS80 +nodefs",
+"+proj=utm +zone=32 +units=m +datum=NAD83 +nodefs",
+"+proj=etmerc +lat_0=0 +y_0=0 +lon_0=9 +x_0=500000 +k=0.9996 +units=m +datum=WGS84",
+"+proj=utm +zone=32 +units=m +ellps=intl +towgs84=-87,-98,-121,0,0,0",
+"+proj=etmerc +lat_0=0 +y_0=-5000000 +lon_0=9 +x_0=200000 +k=0.999980 +units=m +ellps=GRS80"
+]
+
+def test():
+	nerr=0
+	for text in TESTS:
+		mlb=TrLib.FromProj4(text)
+		if mlb is None:
+			mlb="Unable to translate..."
+			ner+=1
+		print("\n%s\n-> %s"%(text,mlb))
+	return nerr
+
+def main(args):
+	progname=os.path.basename(args[0])
+	IS_INIT=TrLib.IS_INIT
+	if not IS_INIT:
+		if "-lib" in args: #In this case we assume that input is the library that we want to test....
+			libpath=args[args.index("-lib")+1]
+			lib=os.path.basename(libpath)
+			dir=os.path.dirname(libpath)
+			IS_INIT=TrLib.InitLibrary("",lib,dir)
+		else:
+			print("You can specify the TrLib-library to use by %s -lib <lib_path>" %progname)
+			IS_INIT=TrLib.InitLibrary("")
+	if not IS_INIT:
+		print("Could not initialize library...")
+		print("Find a proper shared library and geoid dir and try again!")
+		return -1
+	print("Running %s at %s" %(progname,time.asctime()))
+	print("Using TrLib v. '%s'" %TrLib.GetVersion())
+	print("Shared library is: %s" %repr(TrLib.tr_lib))
+	nerr=test()
+	print("Errors: %d" %nerr)
+	return nerr
+
+if __name__=="__main__":
+	sys.exit(main(sys.argv))
+		
+		
  */
  """
 #########################
-## Python bindings for simple test api to KmsTrLib/KMSTrans
-## simlk, may->sep 2011                
+## Python bindings for (Kms) TrLib
+## simlk, 2011,2012            
 #########################
 try:
 	import numpy as np
 import sys
 import time
 IS_INIT=False
-if "win" in sys.platform:
+if sys.platform.startswith("win"):
 	STD_LIB="KMSTRLIB.dll"
+elif "darwin" in sys.platform:
+	STD_LIB="KMSTRLIB.dylib"
 else:
 	STD_LIB="KMSTRLIB.so"
 STD_DIRNAME=os.path.dirname(__file__)

testdata/new/DK_3d_geoHed50_h_dvr90.txt

 9.97661975 55.41850494 36.39012918
 10.20236548 55.37343101 75.30307115
 10.49660915 55.15485477 25.47907175
-10.64183341 55.39415703 82.70625043
+10.64183341 55.39415703 82.70625044
 10.67882413 55.19554863 52.57326743
 10.04956849 55.31455064 18.12353167
 11.14249779 55.69650795 53.97810765
 10.05463574 56.91193181 24.82248208
 10.01315398 56.54202107 19.66176550
 10.69251150 56.39683459 54.83754406
-9.86719981 56.29341173 17.80303945
+9.86719981 56.29341173 17.80303946
 9.41789361 56.26991152 92.01095876
 8.33378539 56.32049162 64.72255269
 8.76991321 55.94318773 18.97712666