gxgy.to.rectindex=function (gx, gy, gridsize = c(20,20), plotdim = c(1000, 500))

{

   badgxgy = (gx < 0 | gy < 0 | gx >= plotdim[1] | gy >= plotdim[2] |

       is.na(gx) | is.na(gy))

   colno = 1 + floor(gx/gridsize[1])

   rowno = 1 + floor(gy/gridsize[2])

   if (length(badgxgy[badgxgy > 0]))

       colno[badgxgy] = rowno[badgxgy] = NA

   return(rowcol.to.rectindex(rowno, colno, gridsize, plotdim))

}

rowcol.to.rectindex=function (rowno, colno, gridsize = c(20,20), plotdim = c(1000, 500))

{

   badrc = (rowno <= 0 | colno <= 0 | rowno > plotdim[2]/gridsize[2] |

       colno > plotdim[1]/gridsize[1])

   rowno = rowno - 1

   colno = colno - 1

   maxrow = floor(plotdim[2]/gridsize[2])

   index = colno * maxrow + rowno + 1

   if (length(badrc[badrc > 0]))

       index[badrc] = NA

   return(index)

}

no = gxgy.to.rectindex(da$gx, data$gy, gridsize = c(20,20), plotdim = plotdim)