The raster package has a function to do this for you in one line. Use layerize():
# make an example raster
r <- raster(nrow=100, ncol=100)
r <- round(runif(ncell(r),1,4),0)
# create presence/absence rasters, stored in a RasterBrick.
r_dummy <- layerize(r)
if grey_raster <= 28523 does not make sense to arcpy. 28523 is the value you want to check against each pixel of grey_raster, which is where Con function comes handy.
In your case, you should use nested conditions like:
arcpy.env.workspace = ndviws
ndvirasterlist = arcpy.ListRasters()
for rimage in ndvirasterlist:
ndvi_raster = arcpy.sa.Raster(...
This is how I would do it (it requires the Spatial Analyst extension which I think you have):
from arcpy.sa import *
arcpy.CheckOutExtension = "Spatial"
ndvi_raster = Raster("ndvivaluefrompixel")
grey_raster = Raster("grayscalepixelvalue")
output_raster = Con(ndvi_raster >= -0.000005 * grey_raster + 0.314367, 2, 1)