Hi list,
Well, this time I’ve a doubt with mapping generation.
I was already able to read and plot shapefiles, plot point on this map. All
this without any problems, but now I want to do something that I think, only
Golden Software Surfer is capable of.
I would like to plot a kriging result on the existing map (map script
below). Well I looked for on the net, but I was not able to be sure of which
packages I’ll need to perform the kriging (maybe geoR?) and then plot it in
the area map.
Is R able to do that? What packages I need to know to do so?
I just need a starting point to go forward.
I need to generate something like this:
http://i288.photobucket.com/albums/ll162/Godrigos/Kriging.jpg
(but without the top right area, just the different grayscale kriging area).
library(xlsReadWrite)
library(maps)
library(mapdata)
library(maptools)
Pontos<-read.xls('Pontos.xls',sheet=1,rowNames=T)
Prof<-read.xls('Pontos.xls',sheet=2,rowNames=T)
Batimetria<-readShapeSpatial(‘Batimetria.shp')
Municipios<-readShapeSpatial('Municipios.shp')
Recifes<-readShapeSpatial('Recifes.shp')
par(mar=c(4,0,0,0))
map('worldHires','brazil',ylim=c(,),xlim=c(,),type='n')
plot(Batimetria,ylab='',xlab='',border='lightgray',add=T)
plot(Recifes,ylab='',xlab='',lty='dashed',add=T)
plot(Municipios,ylab='',xlab='',col=rgb(229,229,229,max=255),add=T)
axis(1,xaxp=c(,,))
axis(2,yaxp=c(,,))
points(Pontos$long,Pontos$lat,pch=20)
text(Pontos$long,Pontos$lat,rownames(Pontos),pos=1,cex=0.5)
text(Prof$Long,Prof$Lat,rownames(Prof),col='darkgray',cex=0.8,srt=-24)
Ps.: Sorry for hiding the maps coordinates. I had to, a demand of my
coworkers.
Thanks in advance.
___________________________________
MSc. <mailto:r.aluizio@gmail.com> Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laboratório de Micropaleontologia
[[alternative HTML version deleted]]
Rodrigo,
This is an old and quite basic Krig, my data was continuous
measurements in lat, long so binned first.
library(geoR)
counts<-bins2d(long, lat,bin=c(0.1,0.1),plot=FALSE, nlevels=15,
color.palette=heat.colors, xaxs='i', yaxs='i', las=1,
main='')
countsgeo<-as.geodata(counts)
xrange<-seq(range(countsgeo$coords[,1])[1],range(countsgeo$coords[,1])
[2], by=bin[1])
yrange<-seq(range(countsgeo$coords[,2])[1],range(countsgeo$coords[,2])
[2], by=bin[2])
loci<-expand.grid(x=range(counts,y=yrange) #krig grid
kc <- krige.conv(countsgeo, loc=loci,
krige=krige.control(cov.pars=c(1, .25))) #can fine
tune
color.palette<-rainbow #use heat etc.
contour(kc,asp=1, col=col<-
c("#FFFFF",color.palette(length(levels)-1)) , add=TRUE) #makes 0
counts white
You could also use image or other 3d graphics.
Hope that Helps,
Jon Loehrke
Graduate Research Assistant
Department of Fisheries Oceanography
School for Marine Science and Technology
University of Massachusetts
200 Mill Road, Suite 325
Fairhaven, MA 02719
jloehrke@umassd.edu
[[alternative HTML version deleted]]
Hi... As u stated earlier that u r able to plot map for kriging.... As i am new to kriging can you please send me the kriging code for a vector file.. Thanks in Advance.... Rodrigo Aluizio wrote:> > Hi list, > > Well, this time I?ve a doubt with mapping generation. > > I was already able to read and plot shapefiles, plot point on this map. > All > this without any problems, but now I want to do something that I think, > only > Golden Software Surfer is capable of. > > I would like to plot a kriging result on the existing map (map script > below). Well I looked for on the net, but I was not able to be sure of > which > packages I?ll need to perform the kriging (maybe geoR?) and then plot it > in > the area map. > > Is R able to do that? What packages I need to know to do so? > > I just need a starting point to go forward. > > > > I need to generate something like this: > > http://i288.photobucket.com/albums/ll162/Godrigos/Kriging.jpg > > (but without the top right area, just the different grayscale kriging > area). > > > > > library(xlsReadWrite) > > library(maps) > > library(mapdata) > > library(maptools) > > Pontos<-read.xls('Pontos.xls',sheet=1,rowNames=T) > > Prof<-read.xls('Pontos.xls',sheet=2,rowNames=T) > > Batimetria<-readShapeSpatial(?Batimetria.shp') > > Municipios<-readShapeSpatial('Municipios.shp') > > Recifes<-readShapeSpatial('Recifes.shp') > > par(mar=c(4,0,0,0)) > > map('worldHires','brazil',ylim=c(,),xlim=c(,),type='n') > > plot(Batimetria,ylab='',xlab='',border='lightgray',add=T) > > plot(Recifes,ylab='',xlab='',lty='dashed',add=T) > > plot(Municipios,ylab='',xlab='',col=rgb(229,229,229,max=255),add=T) > > axis(1,xaxp=c(,,)) > > axis(2,yaxp=c(,,)) > > points(Pontos$long,Pontos$lat,pch=20) > > text(Pontos$long,Pontos$lat,rownames(Pontos),pos=1,cex=0.5) > > text(Prof$Long,Prof$Lat,rownames(Prof),col='darkgray',cex=0.8,srt=-24) > > > > Ps.: Sorry for hiding the maps coordinates. I had to, a demand of my > coworkers. > > > > Thanks in advance. > > ___________________________________ > MSc. <mailto:r.aluizio at gmail.com> Rodrigo Aluizio > Centro de Estudos do Mar/UFPR > Laborat?rio de Micropaleontologia > > > [[alternative HTML version deleted]] > > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/Plotting-a-kriging-on-a-map-tp20842909p21248673.html Sent from the R help mailing list archive at Nabble.com.
