# Bring in data
library(data.table)
library(raster)
set.seed(123)

# Bring in tifs and subset to your 16 states
tifFiles <- list.files(".", pattern=".TIF$", full.names = TRUE)
tifFiles <- paste0("./Suitability_Prob_v",c(30:39,310:315),".TIF")

# How many points do we want to sample from each state?
# I am not sure what is in your csv file, so I sampled some points from a poisson
sampPoints <- data.table(State=tifFiles, N=rpois(length(tifFiles),6))

# Function to read in raster and sample points
readRasterAndSample <- function(rasFile, count) {

    # Read in raster and deal with 0s and NAs
    ras <- raster(rasFile)
    ras[ras < 0] <- 0

    # Convert to data.table
    rasDT <- as.data.frame(ras, xy=TRUE)
    setDT(rasDT)
    setnames(rasDT, c("x","y","values"))

    # Sample non-zero points according to desired count
    # Use raster probability as sample prob
    outPoints <- rasDT[values > 0, .SD[sample(1:.N, size=min(count,.N), prob=values)]]

    # Adjust points so they are the midpoint of the cellgrid
    outPoints[, x := x + 0.5*res(ras)[1]]
    outPoints[, y := y + 0.5*res(ras)[2]]
    outPoints[, state := basename(rasFile)]

    # Combine into list and return
    outList <- list(ras=ras, pts=outPoints)
    return(outList)
}

# Apply function to all states
rasterList <- lapply(1:nrow(sampPoints), function(x)
    with(sampPoints[x], readRasterAndSample(rasFile=State, count=N)))

# Collect all rasters and merge
rasterAll <- do.call(merge, lapply(rasterList, "[[", "ras"))

# Collect all state points and convert to SpatialPointsDataframe
pointDat <- rbindlist(lapply(rasterList, "[[", "pts"), fill=TRUE)
coordinates(pointDat) <- ~x+y
proj4string(pointDat) <- rasterAll@crs@projargs

# Plot to check
plot(rasterAll)
plot(pointDat, add=TRUE)