coxph diagnostics plot for shape of hazard function?

Similarly, when I do plot(zph), B(t) is fairly non-constant.

This isn't inherently a problem for me. I don't need a hard single number
to characterize the shape of the excess risk. However, I'd like to be
able to say
something qualitative about the shape of the excess risk for the predictor.
E.g., is it linear, monotonically increasing, monotonially decreasing, etc.
Is it safe to use the coxph diagnostic plot for this purpose?

Basically - yes you can.  There are a few caveats:
    1. As a computational shortcut cox.zph assumes that var(X) is approximately 
constant over time, where X is the matrix of covariates.  (Improving this has 
been on my to do list for some time). I have found this to be almost always 
true, but if you have a data set where e.g. everyone in treatment 1 is crossed 
over at 6 months, then you can get odd results for that covariate.  I've run 
across 2-3 such data sets in 10+ years.

    2. The spline curve on the plot is "for the eye".  You can certainly use 
other smoothings, fit a line, etc.  Often you can find a simpler fit.
    	zpfit <- cox.zph(mycoxfit, transform='identity')
    	plot(zpfit$x, zpfit$y[,1], xlab='Time')  #look at variable 1
    	lines(lowess(zpfit$x, zpfit$y[,1]), col=2)
    	abline( lm(zpfit$y[,1] ~zpfit$x), col=3)
    	
    	plot(zpfit$x, zpfit$y[,1], log='x') #same as transform=log 
    	
    	etc.
    	
Sometimes the regression spline fit, the default for cox.zph, puts an extra 
"hook" on the end of the curve, somewhat like polynomials will.  
    	
    	Terry T.