R/rawCrop2nicenc.R

In sidjai/biosplit: Model for continental seasonal migration of insects

#'Manipulate the input slices into a nice netCDF file
#'
#'Scales the color tag to a single true or false value,
#'	aggregates the small spacing into the desired spacing using sum(extract()) 
#'	since raster::aggregate can't take uneven spaces.
#'	Then it changes the projection into unprojected lat/lon or a user defined projection 
#'
#'
#'@param dirIn The directory with slices of tif files ussually generated by NASS2TIFs()
#'@param pathOut The full file name of the aggregated file, either in .dat or .nc
#'@param cropGrid The input extent and projection with the new spacing that you want
#'@param niceGrid The result grid in the projection that you want
#'
#'@return saves the result as a csv readable by surfer (.dat) and the netCDF file at pathOut
#'	as specified. Returns a 0 if the writing failed, 1 if it worked.
#'@import ncdf
#'@import raster
#'@export
rawCrop2nicenc <- function(dirIn, pathOut, cropGrid, niceGrid=""){
	
	#Raster definitions start at the top left so need to reverse the list
	tifFiles <- parseFiles(dirIn,'tif',shReverse = TRUE)
	
	xs <- xFromCol(cropGrid)
	aggxs <- list(left = xs - xres(cropGrid)/2,
								right = xs + xres(cropGrid)/2)
	if(is.character(niceGrid)){
		#default out grid is unprojected lat lon
		niceProj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
		niceGrid <- projectExtent(cropGrid, niceProj)
	}
	
	#tests
	
	numSlicesIn <- length(tifFiles)
	if (numSlicesIn != nrow(cropGrid)){
		stop(paste("Definition expects", nrow(cropGrid),
							 "while the tifs in", dirIn, "has", numSlicesIn))
	}
	
	prog <- txtProgressBar(style = 3)
	zs <- vapply(tifFiles,function(x){
		setTxtProgressBar(prog, getTxtProgressBar(prog) + 1/numSlicesIn)
		processCrop(x,aggxs)
	},rep(1,ncol(cropGrid)),USE.NAMES = FALSE)
	
	zs <- as.vector(zs)
	translate <- rasterToPoints(cropGrid)
	
	xyz <- cbind(translate,zs)
	
	
	write.table(xyz,
							file = gsub("[.]nc", ".dat", pathOut),
							quote=FALSE,row.names=FALSE,col.names=FALSE)
	
	ras <- rasterFromXYZ(xyz)
	projection(ras) <- projection(cropGrid)
	ras <- projectRaster(ras,niceGrid)
	dataType(ras) <- 'INT4S'
	sult <- try(writeRaster(ras,filename=gsub("[.]dat", ".nc", pathOut),
													format="CDF", overwrite=TRUE),silent=TRUE)
	
}



processCrop <- function(gdalSlice,defxs){
	ras <- raster(rgdal::readGDAL(gdalSlice, silent = TRUE))
	top <- ras@extent@ymax #get the top of the slice
	bot <- ras@extent@ymin #get the bottom of the slice
	
	#Go through each block and sum them up individually to aggregate
	#raster::aggregate can't do uneven spaces 30 -> 40000
	#0.09 is the conversion from pixel to hectares
	zvals <- vapply(1:length(defxs$left), function(xi){
		ext <- extent(defxs$left[xi], defxs$right[xi], bot, top)
		
		return(0.09*sum(as.logical(extract(ras, ext)),na.rm=TRUE))
	},1)
	ras <- 999
	return(zvals)
}