Water depth - echo sounder (Fahrentholz)
Jump to navigation
Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
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) }