how can I improve my plot and basemap overlay?

Question

I have a basic plot of rainfall contours with a map overlay.
1. I wonder if there is something I can change on the aspect ratio so that my map is proportional. Furthermore I am wondering why the frame disregards the map limits that I set.
2. Is there an extrapolation function I can add so that I can improve my plot and fill the map area with contours, even if it means filling the whole plotting frame.
Help will be appreciated. My code is below.

Libraries

library(maps)  
library(akima)  

Set working directory and get data

setwd("C:/R_Progs/Scripts")  
datr <- read.table("./RFMEAN.DAT",header=TRUE,sep="")  
datfr <- data.frame(datr)  
x <- datfr[,1]  
y <- datfr[,2]  
z <- datfr[,3]  

Interpolation

fld <- with(datfr, interp(x , y, z))  

Make the plot

filled.contour(x = fld$x, y = fld$y, z = fld$z,  
plot.axes = {axis(1); axis(2); map('world', xlim = c(21, 34), ylim =  
c(-19, -8), add = T, col = "darkgrey")},  
color.palette = colorRampPalette(c("white", "blue")),  
xlab = "Longitude",  
ylab = "Latitude",  
main = "Zambian rainfall",  
key.title = title(main = "Rain (mm)", cex.main = 1))  

My data

Lon      Lat   Rain  
32.58   -13.55  1016.9  
27.07   -16.85  801.6  
32.67   -10.17  1090  
24.2    -13.6   1034.4  
28.47   -14.45  907.7  
28.5    -14.4   917.7  
28.12   -12.6   1309.1  
27.92   -15.77  761.4  
22.7    -14.    825.2  
24.8    -14.8   902.4  
31.13   -10.22  1343.2  
25.85   -13.53  1146.8  
29.08    -9.8   1378.6  
25.82   -17.82  691.7  
33.2    -12.28  922.9  
28.32   -15.42  1078.1  
28.45   -15.3   882.1  
27.63   -16.13  737.3  
28.85   -11.1   1176.6  
31.33    -8.85  1267.8  
31.93   -13.27  817.5  
31.22   -10.18  1349.7  
23.15   -15.25  918.2  
31.43   -11.9   1040  
32.56   -13.65  1011.7  
28.25   -15.55  856.7  
27.07   -14.98  891.2  
24.43   -11.75  1413.1  
28.65   -13.    1232.8  
31.28   -14.25  966.8  
29.53   -11.35  1476.5  
23.27   -16.12  729.7  
30.22   -13.23  1133.6  
24.3    -17.47  756.3  
26.38   -12.18  1309.7  
23.12   -13.53  1057.4  

And the plot I get (I hope it is uploaded)

Answer 1

You can't extrapolate outside your data range, so no there is no way to fill the entire map unless you want to make up data.

I like ggmap much better. Maybe you can try/start with this and mess around.

library(ggmap)
z = get_map(location = "Zambia", color = "bw", zoom = 6)
ggmap(z, extent = "normal", maprange=FALSE) +
  stat_density2d(data = df, aes(x = Lon, y = Lat, z = Rain, fill = ..level..),
                 geom="polygon", bins = 100, size = 0.01, alpha = 0.05) +
                 scale_fill_gradient(name = "Rainfall", low = "dodgerblue", high = "dodgerblue4")

Edit for Comment

Run this in case your df is different.

df=structure(list(Lon = c(32.58, 27.07, 32.67, 24.2, 28.47, 28.5, 28.12, 27.92, 22.7, 24.8, 31.13, 25.85, 29.08, 25.82, 33.2, 28.32, 28.45, 27.63, 28.85, 31.33, 31.93, 31.22, 23.15, 31.43, 32.56, 28.25, 27.07, 24.43, 28.65, 31.28, 29.53, 23.27, 30.22, 24.3, 26.38, 23.12), Lat = c(-13.55, -16.85, -10.17, -13.6, -14.45, -14.4, -12.6, -15.77, -14, -14.8, -10.22, -13.53, -9.8, -17.82, -12.28, -15.42, -15.3, -16.13, -11.1, -8.85, -13.27, -10.18, -15.25, -11.9, -13.65, -15.55, -14.98, -11.75, -13, -14.25, -11.35, -16.12, -13.23, -17.47, -12.18, -13.53), Rain = c(1016.9, 801.6, 1090, 1034.4, 907.7, 917.7, 1309.1, 761.4, 825.2, 902.4, 1343.2, 1146.8, 1378.6, 691.7, 922.9, 1078.1, 882.1, 737.3, 1176.6, 1267.8, 817.5, 1349.7, 918.2, 1040, 1011.7, 856.7, 891.2, 1413.1, 1232.8, 966.8, 1476.5, 729.7, 1133.6, 756.3, 1309.7, 1057.4)), .Names = c("Lon", "Lat", "Rain"), class = "data.frame", row.names = c(NA, -36L))