User Tools

Site Tools


teaching:m_r:excersize:proj4

Co-ordinate conversion using Proj4 and Awk

Proj.4

Reference: https://trac.osgeo.org/proj/

PROJ.4 is a library for performing conversions between cartographic projections. The library is based on the work of Gerald Evenden at the USGS,[2] but is now an OSGeo project maintained by Frank Warmerdam. The library also ships with executables for performing these transformations from the command line. (Reference: http://en.wikipedia.org/wiki/PROJ.4)

Download program
WindowsMacLinux
Proj.4 Proj.4

Awk

Awk is a terminal program / script program.

Reference: http://www.gnu.org/software/gawk/

If you are like many computer users, you would frequently like to make changes in various text files wherever certain patterns appear, or extract data from parts of certain lines while discarding the rest. To write a program to do this in a language such as C or Pascal is a time-consuming inconvenience that may take many lines of code. The job is easy with awk, especially the GNU implementation: gawk.

The awk utility interprets a special-purpose programming language that makes it possible to handle simple data-reformatting jobs with just a few lines of code.

Find the manual here: http://www.gnu.org/software/gawk/manual/

Download program
WindowsMacLinux
Gnu Awk Gnu Awk
open gawk.pkg
Gnu Awk

Awk/proj4 script

Script that convert gga-strings to UTM co-ordinates.

ggacs2cs.awk
#!/usr/bin/awk -f
function usage(){
  print "purpose: convert geographical/projected coordinates from NMEA input"
  print "default: parses GGA and transforms coordinates from WGS84 to UTM 32 N"
  print "options [-v OPT=value]:"
  print "EPSGFROM=4326 (epsg code to convert from)"
  print "EPSGTO=25832 (epsg code to convert to)"
  print "MYFS=\",\" (field separator)"
  print "CNV=false (convert to pseudo NMEA message with converted coordinates)"
  print "NFN=3 (input field latitude)"
  print "EFN=5 (input field longitude)"
  print "ECHO=false (echo NMEA messages)"
  print "PRINT2FILE=false (set to false in order not to write to file)"
  print "OUTPUTFILE=ggacs2cs (filename of the output file)"
 
}
#Windows command: gawk -f ggacs2cs.awk -v ECHO=true PRINT2FILE=true OUTPUTFILE=utm_row4_adana EPSGTO=32636 ggaLog.log
#Windows command: gawk -f ggacs2cs.awk -v PRINT2FILE=true OUTPUTFILE=gga_l_u EPSGTO=32636 GGAlinksoben.log
#EPSG:32636: WGS 84 / UTM zone 36N
 
# TODO: feed all coordinates to cs2cs in one large chunk, not line per line
# TODO: use optarg: http://snap.nlc.dcccd.edu/reference/awkref/gawk_17.html
BEGIN{
<<<<<<< .mine
CONVFMT="%.12f"
=======
CONVFMT="%.10f"
>>>>>>> .r1275
#     while ((c = getopt(ARGC, ARGV, "h")) != -1) {
#         if (c == "h") {
#             usage()
#         }# else if (c == "c") { # TODO implement other options here
        #}
 
if(MYFS==""){
  FS=","
}
  OFS=FS
  if (EPSGTO==""){
    EPSGTO="25832" # 25832: # ETRS89 / UTM zone 32N <25832> +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs  <>
  }
  if (EPSGFROM==""){
    EPSGFROM="4326"
  }
  if (CNV==""){ # if CNV="true" then a fake string is created with the transformed coordinates
    CNV="false" 
  }
  if(NFN==""){NFN=3}
  if(EFN==""){EFN=5}
  if(HFN==""){HFN=10}
  if(ECHO==""){ECHO="false"} # echo NMEA as seen? set with -v ECHO=true
  if(PRINT2FILE==""){PRINT2FILE="false"}
  if(OUTPUTFILE==""){OUTPUTFILE="ggacs2cs"}
}
 
function nmea2dd(nmeacoord){
  dn=int(nmeacoord/100);
  return dn+((nmeacoord-(dn*100.0))/60)
}
function cs2cs(coordarray){
  CMD="echo "coordarray[2]" "coordarray[1]" "coordarray[3]" | cs2cs -f \"%.6f\" +init=epsg:"coordarray[4]" +to +init=epsg:"coordarray[5]
   #print CMD
  gkarray[0]=""
  CMD | getline TRANSFCOORD
#   printf "begin"TRANSFCOORD"end"
  split(TRANSFCOORD, gkarray, "\t")
  split(gkarray[2],D," ")
  # print  C[0]
  # GKR=C[1]
  # GKH=D[1]
  coordarray[7]=gkarray[1]
  coordarray[6] = D[1] #northing
  coordarray[8] = D[2] #easting
#   print "\\debug wgs2gk3: " coordarray[1] FS coordarray[2] FS coordarray[3]  FS coordarray[4] FS coordarray[5] FS coordarray[6]
#   return gkarray
}
 
# function wgs84toepsg(coord_n_e_h_to){
#   CMD="echo "coord_n_e_h_to[1]" "coord_n_e_h_to[2]" "coord_n_e_h_to[3]" | cs2cs -f \"%.3f\" +proj=latlong +datum=WGS84 +to +EPSG:"coord_n_e_h_to[4]" "
#   # print CMD
#   resarray[0]=""
#   CMD | getline TRANSFCOORD
# #   printf "begin"TRANSFCOORD"end"
#   split(TRANSFCOORD, splitarray, "\t")
#   split(splitarray[2],D," ")
#   # print  C[0]
#   # GKR=C[1]
#   # GKH=D[1]
#   resarray[1]=splitarray[1]
#   resarray[2] = D[1]
#   resarray[3] = D[2]
#    print "\\debug wgs84toepsg: " resarray[1] FS resarray[2] FS resarray[3]
#   return resarray
# }
{if (ECHO!="false"){print;}
#if (true){
#	for (i=1;i<=NF;i++){
#		print i"\t"$i
#	}
#}
}
/GGA/{
#   $GPGGA,131329.00,4842.71443089,N,00912.27614062,E,4,07,1.3,437.834,M,0.000,M,0.0,0000*70
  transarray[0]=""
  transarray[1]=nmea2dd($NFN) # ddlat
  transarray[2]=nmea2dd($EFN) # ddlon
  transarray[3]=$HFN		# wgsheight
  transarray[4]=EPSGFROM
  transarray[5]=EPSGTO
#   transarray[4]=EPSGCODE
  cs2cs(transarray)
  if(CNV!="false"){ # print pseudo string with transformed coordinates
    $1="$CCONV"
    $NFN=transarray[6]
    $EFN=transarray[7]
    $HFN=transarray[8]
    print $0
  }else
  print transarray[6] FS transarray[7] FS transarray[8]
  #write to file: works in windows
  if(PRINT2FILE=="true"){
	#{print transarray[6] FS transarray[7] FS transarray[8] > OUTPUTFILE".log"}
	#northing,easting
	{print transarray[1] FS transarray[2] FS transarray[3] > OUTPUTFILE".log"}
 
  }
}
teaching/m_r/excersize/proj4.txt · Last modified: 2012/11/07 18:03 by karina