Water depth - echo sounder (Fahrentholz)
Jump to navigation
Jump to search
Parameter to be measured:
Water depth (subsurface measurement)
Method:
echo sounder (sonar)
Equipment:
Fahrentholz BBES 30/ 100 /200 /700
Advantages:
- no-frills sturdy device
- simple recording with any COM-port monitor (e.g. Telnet, Hyperterminal)
- up to 40 measurements/sec from 0.1 to 100 m
- custom modifications on request
- can be used in conjunction with (d)GPS, synchronize both records via timestamp (see scripts below)
- rustic bulky Landrover-style
Disadvantages:
- computer connection needed, no internal logging
- COM-port needed
- external 12/24 V battery needed
- no logging software provided (see scripts below)
- no real-time display of values
What to watch out for:
- may sometimes loose calibration and require resetting (press button and the device is back in seconds)
- when attached to the stern, bubbles in the wake may produce noisy signal
- when used with (d)GPS, synchronize system time of laptop with GPS time. BEWARE: true GPS-time is 15 sec ahead of UTC (Nov 2010), which increases every 8 month due to leap seconds. Check which one the GPS displays and which one it records (this may differ!).
Problems/Questions:
Price:
~ 2500 € (2003, Germany)
Links
Projects that used the above equipment:
Other related web sites:
[[ http://www.fahrentholz.de/bbesneu.htm ]]
Python script (logging with timestamp from serial port, aggregation of values) File:Echolot2.zip
References
Scripts
readme.txt:
echolot ------- reads data from echo-sounder (Fahrentholz Echolot-Sensor) via Serial port, aggregates to specified timestep and writes data and timestamp to file Till Francke; 10.3.2010 Please don't bother me with bugs, improve the script if you can. Possibly better alternative (not tested): http://www.geo.uib.no/polarhovercraft/index.php?n=Main.GPSEchoSounderLogging requirements: - Python (>=2.6.4) with packages numpy (>=1.3.0) and pyserial (>=2.5) - laptop with serial port (or a reliable USB-COM-converter) initial test: - connect to the echo-sounder using a terminal programm (e.g. Hyperterminal, Putty) choosing the right COM-port and 4800 baud settings: - open "echolot.py" with a text editor and customize the following lines. Don't add or remove Tabs. time_interval=3 #averaging interval used for file output [seconds] log_file="e:/temp/comlog.txt" #output file file_mode="w" #file mode: "w": overwrite existing files; "a": append to existing files [...] ser.port = 1 #COM-port 0:com1; 1:com2, etc. - if you want to synchronize the logfile with a GPS-track, verify that the system time is set correctly operation: - run "echolot.py". A console window will open and display current measurements, min and max values - the specified log file should grow postprocessing: - you can use the R-script (needs R installed) "link2GPX.r" as a start for syncing a echo-sounder log with a GPS-track: - convert GPS-log to GPX - adjust filenames near "gps_log=..." and "echo_log=..." - adjust time offset near "time_offsets=c(15,0)". Beware: some GPS-receivers may log a different time than they display. The offset is than the real-time minus GPS-time, which is 15 sec in 2010.
Copy the contents of this box and save as echosound.py. echosound.py:
#!/usr/bin/python #reads data from echo-sounder (Fahrentholz Echolot-Sensor) via Serial port (Com-Port), aggregates and writes to file #Till Francke; 10.3.2010 #needs PySerial, NumPy import serial if 1: time_interval=3 #averaging interval used for file output [seconds] log_file="./comlog.txt" #output file file_mode="a" #file mode: "w": overwrite; "a": append ser = serial.Serial() #ser.baudrate = 57600 #MRR ser.baudrate = 4800 #Echolot ser.port = 0 #COM-port 0:com1; 1:com2 ser.timeout = 0.5 # ser def convertStr(s): try: ret = float(s) except ValueError: ret = NaN return ret import random from datetime import datetime from numpy import * import csv with open( log_file, file_mode) as f: spamWriter = csv.writer(f, delimiter='\t',quotechar='|', lineterminator="\n",quoting=csv.QUOTE_MINIMAL) if file_mode=='w': spamWriter.writerow(["Time","meandepth[cm]","min","max","std"]) #write header ser.open() while 1: now = datetime.today() elapsed = datetime.today()-now errorline="" a=[] #b = array([], dtype=float) while elapsed.seconds < time_interval: line = ser.readline() # read a '\n' terminated line line = line.rstrip() #print("read from port:",line) convertvalue=convertStr(line) if isnan(convertvalue): errorline=line #save uninterpretable line a=a + [convertvalue] # a=a + [random.randint(1, 10)] elapsed = datetime.today()-now print(time_interval," seconds interval") if errorline=="": print 'mean:{0:2f} min:{1:3f} max:{2:4f} std:{3:4f}'.format(mean(a),min(a),max(a),std(a)) spamWriter.writerow([datetime.today().strftime("%Y-%m-%d %H:%M:%S"),mean(a),min(a),max(a),std(a)]) if isnan(mean(a)): print("\a") # beep f.flush() else: print ("received:",errorline, "\a") ser.close() # close port f.close()
Scripts for linking the recodred file to a GPX-Track: link2GPX.r:
##join logger files by timestamp source("join_logfiles_functions.r") #declare necessary functions #import GPS-log gps_log=read.table(file="DGPS_Data_Bathymetrie.csv", sep=";", header=TRUE, dec=",") gps_log$Control=NULL gps_log$Method=NULL gps_log$Orbit=NULL gps_log$Projection=NULL gps_log$time=as.POSIXct(strptime(as.character(gps_log$Start.Time),"%Y-%m-%d %H:%M:%S")) #,tz="GMT" gps_log$Start.Time=NULL gps_log$id=1:nrow(gps_log) #import echo-log echo_log=read.table(file="comlog.txt", sep="\t", header=TRUE) echo_log$time=as.POSIXct(strptime(as.character(echo_log$Time),"%Y-%m-%d %H:%M:%S")) #,tz="GMT" echo_log$Time=NULL ##adjust shift time_diff= diff(as.numeric(gps_log$time)) #compute time difference between succesive log points # time_diffi=which(time_diff==0) # write.table(gps_log[sort(c(time_diffi,time_diffi+1,time_diffi+2,time_diffi-1)),],"duplicates.csv",row.names=F,quote=F,sep="\t") time_diff[time_diff > 5]=0 #discard time differences larger than 5 s gps_log$speed_ms=c(Inf, sqrt(diff(gps_log$Grid.Northing..m.)^2 + diff(gps_log$Grid.Easting..m.)^2) / time_diff) *3.6 gps_log$speed_ms[is.na(gps_log$speed_ms)]=Inf gps_log$speed_ms2=c(gps_log$speed_ms[-1],Inf) #calculate speed from previous and to next GPS-point ##clean gps-log gps_log$Elevation..m.[gps_log$Vertical.Precision..m.>0.5]=Inf #remove elevation records with poor reception # # graphics.off() # firstdaystart=ISOdatetime(2010,10,3,16,17,40) #, tz="GMT" #day1 zoom : 14 s # firstdayend=ISOdatetime(2010,10,3,16,18,30) #, tz="GMT" # # firstdaystart=ISOdatetime(2010,10,7,12,32,15) #, tz="GMT" #day2 zoom : 15 s # firstdayend=ISOdatetime(2010,10,7,12,32,50) #, tz="GMT" # # firstdaystart=ISOdatetime(2010,10,7,15,32,0) #, tz="GMT" #day2 langsam # firstdayend= ISOdatetime(2010,10,7,15,39,0) #, tz="GMT" # # firstdaystart=ISOdatetime(2010,10,7,17,33,40) #, tz="GMT" #day2 ende # firstdayend= ISOdatetime(2010,10,7,17,34,40) #, tz="GMT" # # firstdaystart=ISOdatetime(2010,10,3,16,36,0) #, tz="GMT" #day1 langsam # firstdayend= ISOdatetime(2010,10,3,16,46,0) #, tz="GMT" # # firstdaystart=ISOdatetime(2010,10,3,18,16,0) #, tz="GMT" # day1 ende # firstdayend= ISOdatetime(2010,10,3,20,46,0) #, tz="GMT" # # # first_day_echo=echo_log$time < firstdayend & echo_log$time > firstdaystart # first_day_gps=gps_log$time < firstdayend & gps_log$time > firstdaystart # # plot(echo_log$time[first_day_echo],echo_log$meandepth.cm[first_day_echo]/100,type="l", pch=".", ylab="depth [m] / speed [km/h]", main="day1", xaxt="n", # ylim=c(0, max(echo_log$meandepth.cm[first_day_echo]/100,gps_log$speed_ms[first_day_gps],na.rm=TRUE))) # # axis.POSIXct(1, x=echo_log$time[first_day_echo], format="%H:%M:%S") # lines(gps_log$time[first_day_gps],gps_log$speed_ms[first_day_gps],type="l", pch=".", col="red") # legend("topright",c("depth","speed"), col=c("black","red"),lty=1) # ## echo_log$time <- echo_log$time + 15 # # #plot map # windows() # plot(gps_log$Grid.Easting..m.,gps_log$Grid.Northing..m.,type="p",pch=".") # points(gps_log$Grid.Easting..m.[first_day_gps],gps_log$Grid.Northing..m.[first_day_gps],type="p",pch="o",col="red") # legend("topright",c("all","selection"), col=c("black","red"),lty=1) #join datasets dataset_list=list(echo=echo_log,gps=gps_log) #specify, which dataframes will be joined (the first in list determines the temporal resolution) time_offsets=c(15,0) #this numbers will be added to the respective time_columns to adjust for offsets for (i in 1:length(dataset_list)) dataset_list[[i]]$time <- dataset_list[[i]]$time + time_offsets[i] join_field="time" #field, by which the joining takes place (usually time) options(warn=1) joined_data_set=join_datasets(dataset_list=dataset_list,join_field=join_field) #join the data joined_data_set$gps.speed_ms[joined_data_set$gps.speed_ms==Inf]=NA #this prevents nodata values from being interpolated joined_data_set$gps.speed_ms2[joined_data_set$gps.speed_ms2==Inf]=NA joined_data_set$gps.Elevation..m.[joined_data_set$gps.Elevation..m.==Inf]=NA # analyse dGPS height, correct for water level and speed effects joined_data_set$speed_kmh=apply(joined_data_set[,c("gps.speed_ms","gps.speed_ms2")],1,mean,na.rm=T) campaigns=list( c(start=ISOdatetime(2010,10,3,8,00,00, tz="GMT"), end=ISOdatetime(2010,10,3,19,00,00)), c(start=ISOdatetime(2010,10,7,8,00,00), end=ISOdatetime(2010,10,7,17,34,00)) ) dist_antenna_sensor= 1.635 #distance between echo sensor and GPS antenna in m joined_data_set$bathy_elevation=0 #prepare column for bathymetry information for (i in 1:length(campaigns)) { index=(joined_data_set$time >= campaigns[[i]]["start"] & joined_data_set$time < campaigns[[i]]["end"]) antenna_elev_time=lm(gps.Elevation..m.~time,data=joined_data_set[index,]) antenna_elev_speed=lm(gps.Elevation..m. ~ I(speed_kmh*speed_kmh)+speed_kmh,data=joined_data_set[index,]) antenna_elev_speedtime=lm(gps.Elevation..m. ~ (speed_kmh*time)+I(speed_kmh*speed_kmh),data=joined_data_set[index,],na.action=na.exclude) #time dependence plot(joined_data_set$time[index],joined_data_set$gps.Elevation..m.[index],type="l") timespan=seq(from=min(joined_data_set$time[index]),to=max(joined_data_set$time[index]),length.out=20) lines(timespan, predict(antenna_elev_time,newdata=data.frame(time=timespan)),col="red") antenna_elev_predicted= predict(antenna_elev_speedtime,newdata=data.frame(speed_kmh=mean(joined_data_set$speed_kmh[index],na.rm=T),time=joined_data_set$time[index])) lines(joined_data_set$time[index], antenna_elev_predicted,col="blue") windows() #speed dependence plot(joined_data_set$speed_kmh[index],joined_data_set$gps.Elevation..m.[index],type="p") speedspan=seq(from=min(joined_data_set$speed_kmh[index],na.rm=TRUE),to=max(joined_data_set$speed_kmh[index],na.rm=TRUE),length.out=20) lines(speedspan, predict(antenna_elev_speed,newdata=data.frame(speed_kmh=speedspan)),col="red") antenna_elev_predicted= predict(antenna_elev_speedtime,newdata=data.frame(time=mean(joined_data_set$time[index],na.rm=T),speed_kmh=joined_data_set$speed_kmh[index])) points(joined_data_set$speed_kmh[index], antenna_elev_predicted,col="blue",pch=".") windows() #bivariate dependence antenna_elev_predicted= predict(antenna_elev_speedtime) #plot residuals plot(joined_data_set$time[index],joined_data_set$gps.Elevation..m.[index]-antenna_elev_predicted,type="l") joined_data_set$bathy_elevation[index]=antenna_elev_predicted - dist_antenna_sensor - joined_data_set$meandepth.cm.[index]/100 #correct original gps_log, too index=(gps_log$time >= campaigns[[i]]["start"] & gps_log$time < campaigns[[i]]["end"]) antenna_elev_predicted= predict(antenna_elev_speedtime,newdata=data.frame(time=gps_log$time[index], speed_kmh=apply(gps_log[index,c("speed_ms","speed_ms2")],1,mean,na.rm=T) ,na.rm=T)) gps_log$sensor_elev[index]=antenna_elev_predicted - dist_antenna_sensor } #explore data library(fields) n=round(max(joined_data_set$meandepth.cm., na.rm=TRUE)/100) pal=palette(tim.colors(n=n)) # split.screen( rbind(c(0,.8,0,1), c(.8,1,0,1))) # first image sset=1:nrow(joined_data_set) plot(joined_data_set$gps.Grid.Easting..m.[sset],joined_data_set$gps.Grid.Northing..m.[sset], col= pal[1+floor(joined_data_set$meandepth.cm.[sset]/100)], type="p", pch=".",cex=2) if (exists('point_pairs')) for (i in (1:nrow(point_pairs))) { pair=which(joined_data_set$gps.id %in% point_pairs[i,]) if (length(pair)==2) { lines (joined_data_set$gps.Grid.Easting..m.[pair],joined_data_set$gps.Grid.Northing..m.[pair], type="l", col="black") points(joined_data_set$gps.Grid.Easting..m.[pair],joined_data_set$gps.Grid.Northing..m.[pair], type="p", pch="o",cex=2,col="black") } } # screen(2) # plot(rep(1,n),1:n, type='p', pch=22,lwd=2, col=pal, bg=pal, xaxt="n",xlab="", ylab="depth (m)", oma=5+c(1,2,3,4)) # image.plot(matrix(joined_data_set$meandepth.cm.,ncol=1),legend.only=TRUE, smallplot=c(0,.3, .3,.7), # col=pal,main='%') #colorbar image.plot(matrix(joined_data_set$meandepth.cm.,ncol=1),legend.only=TRUE, smallplot=c(0.0,.1, .1,.7), col=pal,main='%') #colorbar # screen(1) # identify(joined_data_set$gps.Grid.Easting..m.[sset],joined_data_set$gps.Grid.Northing..m.[sset],joined_data_set$meandepth.cm.[sset]) close.screen(all=TRUE) source("zoom_plot.r") zoom_plot(joined_data_set$gps.Grid.Easting..m.,joined_data_set$gps.Grid.Northing..m.,joined_data_set$meandepth.cm.) identify(joined_data_set$gps.Grid.Easting..m.,joined_data_set$gps.Grid.Northing..m.,joined_data_set$meandepth.cm.) #find adjacent points search_radius=sqrt(1) #maximum spatial distance in m min_time = 2*60 # minimum temporal distance in s num_time= as.numeric(gps_log$time) #speed point_pairs=NULL for (i in 1:nrow(gps_log)) { valid_points=(gps_log$id > i) & (abs(num_time[i]-num_time)>=min_time) #determine points not yet treated and with minimum temporal distance to current point dists2 = ( gps_log$Grid.Easting..m.[valid_points] - gps_log$Grid.Easting..m.[i] )^2 + ( gps_log$Grid.Northing..m.[valid_points] - gps_log$Grid.Northing..m.[i] )^2 within_dist= dists2<=search_radius if (any(within_dist)) { closest=which.min(dists2) point_pairs=rbind(point_pairs,c(i, which(valid_points)[closest])) print(i); flush.console() } } save(file="pointpairs.Rdata", list="point_pairs") #use replicate points only correlations=NULL gps=gps_log[unique(as.vector(point_pairs)),] for (offsettime in -5:20) { echo=echo_log echo$time=echo$time+offsettime echo= echo[echo$time %in% intersect(echo$time,gps$time),] dataset_list=list(echo=echo,gps=gps) #specify, which dataframes will be joined (the first in list determines the temporal resolution) join_field="time" #field, by which the joining takes place (usually time) options(warn=1) joined_data_set=join_datasets(dataset_list=dataset_list,join_field=join_field) #join the data #plot comparison of depth values of intersection points depths=NULL sd_depth=NULL hor_prec_gps=NULL nomatch=(joined_data_set$gps.id != round(joined_data_set$gps.id) | (abs(joined_data_set$gps.bathy_elevation)==Inf)) joined_data_set=joined_data_set[!nomatch,] for (i in 1:nrow(point_pairs)) { pair_rows=joined_data_set$gps.id %in% point_pairs[i,] # coords=t(joined_data_set[pair_rows,"meandepth.cm."]) coords=t(joined_data_set[pair_rows,c("gps.sensor_elev","meandepth.cm.")]) if (ncol(coords) != 2) next else depths=rbind(depths,coords[1,]-coords[2,]/100) sd_depth=c(sd_depth,sum(joined_data_set[pair_rows,"std"])) hor_prec_gps=c(hor_prec_gps,sum(joined_data_set[pair_rows,"gps.Horizontal.Precision..m."])) } # plot(depths,xlab="depth 1st pass [cm]",ylab="depth 2nd pass [cm]",main="offset 15 s", col= gray((hor_prec_gps-min(hor_prec_gps))/(max(hor_prec_gps)-min(hor_prec_gps))), type="p", pch=".",cex=10) plot(depths,xlab="depth 1st pass [cm]",ylab="depth 2nd pass [cm]",main="offset 15 s, speed corrected", col= "magenta", type="p", pch="o") abline(0,1) correlations=rbind(correlations,c(offsettime, cor(depths[,1],depths[,2],use="pairwise.complete.obs"))) } plot(correlations[,1],correlations[,2], ,xlim=c(-5,20),xlab="temporal offset of echo log [min]",ylab="R2 between depths at same location") #export column_order=c( #for ordering the columns to match for import as 3d-ASCII # which(names(joined_data_set) == "gps.id"), which(names(joined_data_set) == "gps.Grid.Northing..m."), which(names(joined_data_set) == "gps.Grid.Easting..m."), which(names(joined_data_set) == "bathy_elevation") # which(names(joined_data_set) == "time") ) #for import into RealWorks write.table(file="joined_data.asc",na.omit(joined_data_set[,column_order]),quote=FALSE,row.names=FALSE,sep="\t",col.names=F) #write joined data to file #all data write.table(file="joined_data.txt",joined_data_set,quote=FALSE,row.names=FALSE,sep="\t") #write joined data to file
join_logfiles_functions.r:
##auxiliary functions for joining of logfiles (e.g. by timestamp) parse_gpx= function (gpxfile) #extract lat, lon and time from all waypoints and trackpoint in GPX-file { #gpx=scan(file = gpxfile, what = character(), n = -1) #gpx=paste(gpx,collapse=" ",sep=" ") zz <- file(gpxfile, "rb") gpx=readChar(zz, 1e8, useBytes = FALSE) close(zz) extract_value=function(haystack, needle) #extract the value from structure blabla needle=value { keyval=gregexpr(paste(needle,"[^\"]*",sep=""), chunk, ignore.case = TRUE,perl=T)[[1]] #extract key and value val=sub(needle,"",substr(haystack,keyval[1],keyval[1]+attr(keyval,"match.length")-1)) #extract value return(val) } extract_value2=function(haystack, needle) #extract the value from structure blabla <needle> value </needle> { keyval=gregexpr(paste("<",needle,">[^<]*",sep=""), chunk, ignore.case = TRUE,perl=T)[[1]] #extract key and value val=sub(paste("<",needle,">",sep=""),"",substr(haystack,keyval[1],keyval[1]+attr(keyval,"match.length")-1)) #extract value return(val) } pointframe=data.frame(ptype=rep("p",0),lat=rep("p",0),lon=rep("p",0),time=rep("p",0),stringsAsFactors=F) pointpos=gregexpr("<wpt|<trkpt", gpx, ignore.case = TRUE)[[1]] pointposwithend= c(pointpos, nchar(gpx)) for (i in 1:length(pointpos)) { chunk=substr(gpx,pointpos[i],pointposwithend[i+1]-1) #treat substring only if (substr(chunk,1,4)=="<wpt") ptype="wpt" else ptype="tpt" #determine point type pointframe=rbind(pointframe,data.frame(ptype=ptype,lat=extract_value(chunk, "lat[ ]*=[ ]*\""),lon=extract_value(chunk, "lon[ ]*=[ ]*\""),time=extract_value2(chunk, "time[ ]*"),stringsAsFactors=F)) } return(pointframe) } join_datasets=function(dataset_list,join_field) #join numeric columns of datasets contained in the named list "dataset_list" using the column "join_field", which must exist in all datasets #the column "join_field" of the first set is used in the resulting output dataset, all other values are linearly interpolated { for (i in 1:length(dataset_list)) #check availability of time information if (!(join_field %in% names(dataset_list[[i]]))) stop(paste("Dataset",i,"doesn't have a column",join_field)) joined_data_set=dataset_list[[1]] #use first dataset as target resolution for (i in 2:length(dataset_list)) #successively join all datasets { for (j in 1:ncol(dataset_list[[i]])) { if (names(dataset_list[[i]])[j] == join_field) next #skip join field if (is.numeric(dataset_list[[i]][,j])) method="linear" else next #only numbers can be interpolated column = approx(x=dataset_list[[i]][,join_field], y=dataset_list[[i]][,j], joined_data_set[,join_field], rule=c(2,2), method=method) #linear interpolation of new column # browser() joined_data_set=cbind(joined_data_set,column$y) #join interpolated column to common dataset names(joined_data_set)[ncol(joined_data_set)] = paste(names(dataset_list)[i],".",names(dataset_list[[i]])[j],sep="") #label column } } return(joined_data_set) }
plot_functions.r:
collect_legend=function(namestr,pch,lty,col) { assign("namestr_all", cbind(namestr_all,namestr), env = .GlobalEnv) assign("pch_all", cbind(pch_all,pch), env = .GlobalEnv) assign("lty_all", cbind(lty_all,lty), env = .GlobalEnv) assign("col_all", cbind(col_all,col), env = .GlobalEnv) } reset_legend=function() { assign("namestr_all", NULL, env = .GlobalEnv) assign("pch_all", NULL, env = .GlobalEnv) assign("lty_all", NULL, env = .GlobalEnv) assign("col_all", NULL, env = .GlobalEnv) } draw_legend=function(xpos,ypos) { legend(xpos, ypos, namestr_all, col=col_all, pch=pch_all,lty=lty_all) } #break string into substring of maxlenght length using newline chars str_break = function (stri,minlength,maxlength) { sparts=strsplit(stri, " ") sparts=sparts[[1]] #sapply(sparts,nchar) resultstr="" linestr="" while (!is.na(sparts[1])) { len=nchar(sparts[1]) if (len+nchar(linestr)<=maxlength) #substring conveniently fits into new line { linestr=paste(linestr,sparts[1]," ",sep="") #append substring to current line sparts=sparts[-1] #mark this substr as treated } else if (len < maxlength) #substring fits completely into new line { resultstr=paste(resultstr,linestr,"\n",sep="") #append current line to resultstring linestr=sparts[1] #empty current line and append substring to current line sparts=sparts[-1] #mark this substr as treated } else { #substring is too long anyway, use the knife limit=maxlength-nchar(linestr) sparts=c(substr(sparts[1],1,limit),substr(sparts[1],limit+1,1000000),sparts) sparts=sparts[-3] } if (nchar(linestr)>=maxlength-1) #flush the rest { resultstr=paste(resultstr,linestr,"\n",sep="") #append current line to resultstring linestr="" } } resultstr=paste(resultstr,linestr) #append current line to resultstring return (resultstr) } #return date-vector for reasonable scaling and labelling of time-axis #time_range: vector containing min and max of desired timespan dateticks=function(time_range) { mintime=as.numeric(min(time_range))+strptime("1970-01-01", "%Y-%m-%d", tz="GMT") maxtime=as.numeric(max(time_range))+strptime("1970-01-01", "%Y-%m-%d", tz="GMT") time_range=difftime(maxtime,mintime,units="hours") if (time_range<48) { by="1 hours" format="%d/%m %Hh" } else if (time_range<7*24) { by="1 days" format="%d/%m/%y" } else if (time_range<30*24) { by="2 days" format="%d/%m/%y" } else if (time_range<180*24) { by="1 weeks" format="%d/%m/%y" } else if (time_range<366*24) { by="2 weeks" format="%d/%m/%y" } else if (time_range<3*366*24) { by="1 months" format="%m/%y" } else if (time_range<10*366*24) { by="6 months" format="%m/%y" } else { by="1 years" format="%y" } at = seq(mintime, maxtime, by=by) #define tick locations return(list(at,format)) }
zoom_plot.r:
#use zooming/panning functions #2009-04-06 #Till: dynamic labelling of time axis ################################################################################ ################################################################################ ################################################################################ source("plot_functions.r") zoom <- function(plotfunc, extent, xdata,ydata,zdata){ # (C) 2007 GPL by Huidae Cho <http://geni.ath.cx> # Simple R function for interactive zooming/panning # # Example: # data <- runif(100)*10 # extent <- list(x=c(1, 100), y=c(0, 10)) # plotfunc <- function(lim){ # plot(data, xlim=lim$x, ylim=lim$y) # abline(mean(data), 0, col="red") # } # zoom(plotfunc, extent) print("Zoom in: Click two corners", quote=FALSE) print("Zoom out: Click above plot", quote=FALSE) print("Prev extent: Click left of plot", quote=FALSE) print("Next extent: Click right of plot", quote=FALSE) print("Full extent: Click below plot", quote=FALSE) print("Pan: Double click", quote=FALSE) print("Quit: Right button", quote=FALSE) lim <- extent lim.stack <- c(lim$x, lim$y) lim.depth <- 1 lim.cur <- 1 repeat{ print(lim) flush.console() plotfunc(lim,xdata,ydata,zdata) options(locatorBell=FALSE) l <- locator(1) if(is.null(l)) break ext <- par()$usr if(l$x < ext[1] || l$x > ext[2]){ cur <- lim.cur lim.cur <- if(l$x < ext[1]) max(lim.cur-1, 1) else min(lim.cur+1, lim.depth) if(lim.cur != cur) lim <- list(x=lim.stack[lim.cur, 1:2], y=lim.stack[lim.cur, 3:4]) next } if(l$y < ext[3]) lim <- extent else if(l$y > ext[4]){ cx <- (lim$x[1] + lim$x[2]) / 2 cy <- (lim$y[1] + lim$y[2]) / 2 w <- lim$x[2] - lim$x[1] h <- lim$y[2] - lim$y[1] lim <- list(x=c(cx-w, cx+w), y=c(cy-h, cy+h)) }else{ l2 <- locator(1) if(is.null(l2)) break if(sum(l$x == l2$x) || sum(l$y == l2$y)){ w <- lim$x[2] - lim$x[1] h <- lim$y[2] - lim$y[1] lim <- list(x=c(l2$x-w/2, l2$x+w/2), y=c(l2$y-h/2, l2$y+h/2)) }else lim <- list(x=sort(c(l$x, l2$x)), y=sort(c(l$y, l2$y))) } if(lim.cur < lim.depth){ lim.stack <- lim.stack[-((lim.cur+1):lim.depth),] lim.depth <- lim.cur } lim.stack <- rbind(lim.stack, c(lim$x, lim$y)) lim.depth <- lim.depth + 1 lim.cur <- lim.cur + 1 } lim } ################################################################################ ################################################################################ ################################################################################ # The function that plots the time series plotfunc <- function(lim,xdata,ydata,zdata) { scales=max(diff(lim$x), diff(lim$y)) # Plot the series # first image #screen(1) # browser() plot(xdata,ydata, col= pal[round(zdata/100)], type="p", pch=".",cex=max(2,3000/scales), xlim=lim$x, ylim=lim$y) if (exists('point_pairs')) for (i in (1:nrow(point_pairs))) { pair=which(joined_data_set$gps.id %in% point_pairs[i,]) if (length(pair)==2) { lines (joined_data_set$gps.Grid.Easting..m.[pair],joined_data_set$gps.Grid.Northing..m.[pair], type="l", col="black") points(joined_data_set$gps.Grid.Easting..m.[pair],joined_data_set$gps.Grid.Northing..m.[pair], type="p", pch="o",cex=2,col="black") } } } zoom_plot = function(xdata,ydata,zdata) { n=round(max(zdata, na.rm=TRUE)/100) pal=palette(tim.colors(n=n)) # Set full plot extent extent= list(x=range(xdata),y=range(ydata)) # Plot with zoom option # windows(width=7,height=5) split.screen( rbind(c(0,.8,0,1), c(.8,1,0,1))) screen(2) image.plot(matrix(joined_data_set$meandepth.cm.,ncol=1),legend.only=TRUE, smallplot=c(0,.3, .3,.7), col=pal,main='%') #colorbar screen(1) zoom(plotfunc, extent,xdata,ydata,zdata) close.screen(all=TRUE) }