Water depth - echo sounder (Fahrentholz)
(Redirected from Water depth - echo sound (Farentholz))
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)
}