Felipe, boa tarde!

Segue uma sugestão, tentando preservar o máximo do seu código original. Alterações indicadas no script.

Para diminuir o efeito de "serrilhado" nas figuras, aumente a resolução da grade.

### <code r>


# ==== Carregando pacotes necessários ================ 
sapply(c("sp", "maptools", "spdep", "mgcv", "rgdal"), require, char=T)
setwd("D:/Temp/TCC")

##### lendo dados
cwb<-readOGR("div_municipal.shp",layer="div_municipal") 
neo <- read.table('neo_contr1.txt',h=T) 

### Plotando o mapa de Curitiba 
plot(cwb) 
bbox(cwb) 
### Inserindo os controles 
points(neo$XCOORD[neo$Y==0],neo$YCOORD[neo$Y==0],col='darkgreen',cex = 0.5) 
### Inserindo os caso 
points(neo$XCOORD[neo$Y==1],neo$YCOORD[neo$Y==1],col='red',cex = 0.5) 

###### Ajustando modelo GAM 
models <- gam(neo$Y ~ +s(XCOORD,YCOORD, k=10, bs="tp"),family=binomial,data=neo) 
summary(models) 

## Fazendo a grid 
xLim <- range(pretty(bbox(cwb)[1,])) ### Adicionado
yLim <- range(pretty(bbox(cwb)[2,])) ### Adicionado
gx <- seq(xLim[1],  xLim[2], 500); nx <- length(gx); nx ### Modificado
gy <- seq(yLim[1],  yLim[2], 500); ny <- length(gy); ny ### Modificado
gr <- expand.grid(gx, gy) 
XCOORD <- gr$Var1 
YCOORD <- gr$Var2 

### criando data frame para predição 
ndados <- data.frame(XCOORD, YCOORD); head(ndados) 
pred   <- predict(models,ndados,type="response",se.fit=TRUE) 
pred   <- data.frame(pred) 
plot(cwb) 

## image 
image(gx,gy,matrix(pred$fit,ncol=ny,nrow=nx), asp=T, col=heat.colors(3,alpha=.5), ### Modificado
       main="Thin Plate Splines com 128 pontos", 
       xlab="Coordenada X", ylab="Coordenada Y" ) 

plot(cwb, add=T) 
#points(neo$XCOORD,neo$YCOORD,pch=19,cex=0.5,col="white") 

##vis.gam 
vis.gam(models,plot.type="contour",view=c('XCOORD','YCOORD'),asp=1,ylim=c(7160000,7200000), 
main="Usando Thin Plate", 
color='cm') 
plot(cwb,add=T) 

### Adicionado ###############################################################
### Criando uma Máscara - pontos da grade vs. polígono do município
gr.sp   <- gr; coordinates(gr.sp) <- ~Var1+Var2
cwb.pol <- SpatialPolygons(list(Polygons(list(Polygon(cwb@lines[[1]]@Lines[[1]]@coords)), "CWB")))
mask    <- over(gr.sp, cwb.pol)
image(gx, gy, matrix(mask, nrow=nx, ncol=ny), asp=T, axes=T, add=F)

### Usando a máscara
masked <- matrix(pred$fit, nrow=nx, ncol=ny)* matrix(mask, nrow=nx, ncol=ny)

par(mfrow=c(1,2))

image(gx, gy, masked, asp=T, col=heat.colors(3,alpha=.5), ### Modificado
       main="Thin Plate Splines com 128 pontos", 
       xlab="Coordenada X", ylab="Coordenada Y" ); plot(cwb, add=T) 

contour(gx, gy, masked, asp=T, col=3, ### Modificado
       main="Thin Plate Splines com 128 pontos", 
       xlab="Coordenada X", ylab="Coordenada Y" ); plot(cwb, add=T) 
### </code>






================================================
Éder Comunello
Agronomist (UEM), MSc in Environ. Sciences (UEM)
DSc in Agricultural Systems Engineering (USP/Esalq)
Brazilian Agricultural Research Corporation (Embrapa)
Dourados, MS, Brazil |<O>|
================================================
GEO, -22.2752, -54.8182, 408m
UTC-04:00 / DST: UTC-03:00