diff --git a/Makefile b/Makefile
index 23ef504..56926bb 100644
--- a/Makefile
+++ b/Makefile
@@ -8,7 +8,7 @@ LIBS		= -lm -lpthread -ldl
 
 VPATH		= models
 objects 	= main.o cost.o ecc33.o ericsson.o fspl.o hata.o itwom3.0.o \
-		  los.o sui.o pel.o egli.o inputs.o outputs.o image.o image-ppm.o
+		  los.o sui.o pel.o egli.o inputs.o outputs.o image.o image-ppm.o tiles.o
 
 GCC_MAJOR	:= $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 1)
 GCC_MINOR	:= $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 2)
@@ -52,6 +52,8 @@ image-ppm.o: image-ppm.cc
 los.o: los.cc common.h main.hh cost.hh ecc33.hh ericsson.hh fspl.hh hata.hh \
        itwom3.0.hh sui.hh pel.hh egli.hh
 
+tiles.o: tiles.cc tiles.hh common.h
+
 .PHONY: clean
 clean:
 	rm -f $(objects) signalserver signalserverHD signalserverLIDAR
diff --git a/common.h b/common.h
index 7abd1ad..73568a2 100644
--- a/common.h
+++ b/common.h
@@ -22,6 +22,8 @@
 #define	KM_PER_MILE	1.609344
 #define FOUR_THIRDS	1.3333333333333
 
+#define MAX(x,y)((x)>(y)?(x):(y))
+
 struct dem {
 	float min_north;
 	float max_north;
diff --git a/inputs.cc b/inputs.cc
index 1f45304..00ebceb 100644
--- a/inputs.cc
+++ b/inputs.cc
@@ -7,7 +7,7 @@
 #include <limits.h>
 #include "common.h"
 #include "main.hh"
-
+#include "tiles.hh"
 
 int loadClutter(char *filename, double radius, struct site tx)
 {
@@ -125,105 +125,19 @@ int loadClutter(char *filename, double radius, struct site tx)
 	return 0;
 }
 
-void readLIDAR(FILE *fd, int h, int w, int indx,double n, double e, double s, double west)
-{
-	int x = 0, y = 0, reads = 0, a=0, b=0, avg=0, tWidth = 0, tHeight = 0;
-	char line[25000];
-	char *pch;
-
-	dem[indx].max_north=n;
-	dem[indx].min_west=e;
-	dem[indx].min_north=s;
-	dem[indx].max_west=west;
-
-	if (max_north == -90)
-		max_north = dem[indx].max_north;
-
-	else if (dem[indx].max_north > max_north)
-		max_north = dem[indx].max_north;
-
-	if (min_north == 90)
-		min_north = dem[indx].min_north;
-
-	else if (dem[indx].min_north < min_north)
-		min_north = dem[indx].min_north;
-
-	if (dem[indx].max_west > max_west)
-		max_west = dem[indx].max_west;
-	if (dem[indx].min_west < min_west)
-		min_west = dem[indx].min_west;
-
-	if (max_west == -1) {
-		max_west = dem[indx].max_west;
-	} else {
-		if (abs(dem[indx].max_west - max_west) < 180) {
-			if (dem[indx].max_west > max_west)
-				max_west = dem[indx].max_west;
-		} else {
-			if (dem[indx].max_west < max_west)
-				max_west = dem[indx].max_west;
-		}
-	}
-
-	if (min_west == 360) {
-		min_west = dem[indx].min_west;
-	} else {
-		if (fabs(dem[indx].min_west - min_west) < 180.0) {
-			if (dem[indx].min_west < min_west)
-				min_west = dem[indx].min_west;
-		} else {
-			if (dem[indx].min_west > min_west)
-				min_west = dem[indx].min_west;
-		}
-	}
-
-	for (y = h-1; y > -1; y--) {
-		x = w-1;
-	
-		if (fgets(line, 25000, fd) != NULL) {
-				pch = strtok(line, " "); // split line into values
-
-				while (pch != NULL && x > -1) {
-					if (atoi(pch) <= -9999)
-						pch = "0";
-					dem[indx].data[y][x] = atoi(pch);
-					dem[indx].signal[x][y] = 0;
-					dem[indx].mask[x][y] = 0;
-					if (atoi(pch) > dem[indx].max_el) {
-						dem[indx].max_el = atoi(pch);
-						max_elevation = atoi(pch);
-					}
-					if (atoi(pch) < dem[indx].min_el) {
-						dem[indx].min_el = atoi(pch);
-						min_elevation = dem[indx].min_el;
-					}
-
-   					x--;
-					pch = strtok(NULL, " ");
-				}//while
-
-
-		} else {
-			fprintf(stderr, "LIDAR error @ x %d y %d indx %d\n",
-					x, y, indx);
-		}//if
-	}//for
-
-}
-
-int loadLIDAR(char *filenames)
+int loadLIDAR(char *filenames, int resample)
 {
 	char *filename;
 	char *files[100]; // 10x10 tiles
-	int x, y, indx = 0, fc = 0, hoffset = 0, voffset = 0, pos,
-	    dem_alloced = 0;
-	double xll, yll, xur, yur, cellsize, avgCellsize;
+	int indx = 0, fc = 0, hoffset = 0, voffset = 0, pos, success;
+	double xll, yll, xur, yur, cellsize, avgCellsize = 0, smCellsize = 0;
 	char found, free_page = 0, jline[20], lid_file[255],	
-	     path_plus_name[255], *junk = NULL;
+	path_plus_name[255], *junk = NULL;
 	char line[25000];
 	char * pch;
-    	double TO_DEG = (180 / PI);
+    double TO_DEG = (180 / PI);
 	FILE *fd;
+	tile_t *tiles;
 
 	// test for multiple files
 	filename = strtok(filenames, " ,");
@@ -233,113 +147,217 @@ int loadLIDAR(char *filenames)
 		fc++;
 	}
 
-	while (indx < fc) {
-		if( (fd = fopen(files[indx], "rb")) == NULL )
-			return errno;
+	/* Allocate the tile array */
+	if( (tiles = (tile_t*) calloc(fc, sizeof(tile_t))) == NULL )
+		return ENOMEM;
 
-		if (fgets(line, 255, fd) != NULL) {
-			pch = strtok (line," ");
-			pch = strtok (NULL, " ");
-			width = atoi(pch); // ncols
+	/* Load each tile in turn */
+	for (indx = 0; indx < fc; indx++) {
 
-			if (debug) {
-				fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, width);
-				fflush(stderr);
-			}
-
-                    if (fgets(line, 255, fd) != NULL)
-                            height = atoi(pch); // nrows
-
-			if (!dem_alloced) {
-				//Reduce MAXPAGES to increase speed
-				MAXPAGES=fc;
-				if(width>height){
-					IPPD = width;
-				}else{
-					IPPD = height;
-				}
-				// add fudge as reprojected tiles sometimes vary by a pixel or ten
-				IPPD+=50;
-				ARRAYSIZE = (MAXPAGES * IPPD)+50;
-				do_allocs();
-				dem_alloced = 1;
-			}
-
-		}
-		if (fgets(line, 255, fd) != NULL) {
-			sscanf(pch, "%lf", &xll); // xll
+		/* Grab the tile metadata */
+		if( (success = tile_load_lidar(&tiles[indx], files[indx])) != 0 ){
+			fprintf(stderr,"Failed to load LIDAR tile %s\n",files[indx]);
+			fflush(stderr);
+			free(tiles);
+			return success;
 		}
 
-		if (fgets(line, 255, fd) != NULL) {
-			sscanf(pch, "%lf", &yll); // yll
+		if (debug) {
+			fprintf(stderr, "Loading \"%s\" into page %d with width %d...\n", files[indx], indx, tiles[indx].width);
+			fflush(stderr);
 		}
 
-		if (fgets(line, 255, fd) != NULL)
-			sscanf(pch, "%lf", &cellsize); 
+		// Increase the "average" cell size
+		avgCellsize += tiles[indx].cellsize;
+		// Update the smallest cell size
+		if (smCellsize == 0 || tiles[indx].cellsize < smCellsize) {
+			smCellsize = tiles[indx].cellsize;
+		}
 
-		avgCellsize=avgCellsize+cellsize;
+		// Update a bunch of globals
+		if (tiles[indx].max_el > max_elevation)
+			max_elevation = tiles[indx].max_el; 
+		if (tiles[indx].min_el < min_elevation)
+			min_elevation = tiles[indx].min_el;
 
-		/*if(cellsize>=0.5){ // 50cm LIDAR?
-			// compute xur and yur with inverse haversine if cellsize in *metres*
-			double roundDistance = (width*cellsize)/6371000;
-			yur = asin(sin(yll*DEG2RAD) * cos(roundDistance) + cos(yll * DEG2RAD) * sin(roundDistance) * cos(0)) * TO_DEG;
-			xur = ((xll*DEG2RAD) + atan2(sin(90*DEG2RAD) * sin(roundDistance) * cos(yll*DEG2RAD), cos(roundDistance) - sin(yll * DEG2RAD) * sin(yur*DEG2RAD))) * TO_DEG;
-		}else{*/
-			// Degrees with GDAL option: -co "FORCE_CELLSIZE=YES"
-			xur = xll+(cellsize*width);
-			yur = yll+(cellsize*height);
-		//}
+		if (max_north == -90 || tiles[indx].max_north > max_north)
+			max_north = tiles[indx].max_north;
 
-		if (xur > eastoffset)
-			eastoffset = xur;
-		if (xll < westoffset)
-			westoffset = xll;
+		if (min_north == 90 || tiles[indx].min_north < min_north)
+			min_north = tiles[indx].min_north;
 
-		if (debug)
-			fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur);
+		if (tiles[indx].max_west > max_west)
+			max_west = tiles[indx].max_west;
+		if (tiles[indx].min_west < min_west)
+			min_west = tiles[indx].min_west;
 
-
-		// Greenwich straddling hack
-		if (xll <= 0 && xur > 0) {
-			xll = (xur - xll); // full width
-			xur = 0.0; // budge it along so it's west of greenwich
-			delta = eastoffset; // add to Tx longitude later
+		if (max_west == -1) {
+			max_west = tiles[indx].max_west;
 		} else {
-			// Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;)
-			if (xll >= 0)
-				xll = 360-xll;
-			if(xur >= 0)
-				xur = 360-xur;
-			if(xll < 0)
-				xll = xll * -1;
-			if(xur < 0)
-				xur = xur * -1;
+			if (abs(tiles[indx].max_west - max_west) < 180) {
+				if (tiles[indx].max_west > max_west)
+					max_west = tiles[indx].max_west;
+			} else {
+				if (tiles[indx].max_west < max_west)
+					max_west = tiles[indx].max_west;
+			}
 		}
-		if (debug)
-			fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", yll, yur, xur, xll, delta);
+
+		if (min_west == 360) {
+			min_west = tiles[indx].min_west;
+		} else {
+			if (fabs(tiles[indx].min_west - min_west) < 180.0) {
+				if (tiles[indx].min_west < min_west)
+					min_west = tiles[indx].min_west;
+			} else {
+				if (tiles[indx].min_west > min_west)
+					min_west = tiles[indx].min_west;
+			}
+		}
+
+	}
+
+	/* Iterate through all of the tiles to find the smallest resolution. We will
+	 * need to rescale every tile from here on out to this value */
+	int smallest_res = 0;
+	for (size_t i = 0; i < fc; i++) {
+		if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){
+			smallest_res = tiles[i].resolution;
+		}
+	}
+
+	/* Now we need to rescale all tiles the the lowest resolution or the requested resolution. ie if we have
+	 * one 1m lidar and one 2m lidar, resize the 2m to fake 1m */
+	int desired_resolution = resample != 0 && smallest_res < resample ? resample : smallest_res;
+	if (desired_resolution > resample && debug )
+		fprintf(stderr, "Warning: Unable to rescale to requested resolution\n");
+	for (size_t i = 0; i< fc; i++) {
+		float rescale = tiles[i].resolution / (float)desired_resolution;
+		if (rescale != 1){
+			if( (success = tile_rescale(&tiles[i], rescale) != 0 ) ){
+				fprintf(stderr, "Error resampling tiles\n");
+				return success;
+			}
+
+		}
+	}
+
+	/* Now we work out the size of the giant lidar tile. */
+	double total_width = max_west - min_west >= 0 ? max_west - min_west : max_west + (360 - min_west);
+	double total_height = max_north - min_north;
+	if (debug) {
+		fprintf(stderr, "totalwidth: %.7f - %.7f = %.7f\n", max_west, min_west, total_width);
+		fprintf(stderr,"mw:%lf Mnw:%lf\n", max_west, min_west);
+	}
+	/* This is how we should _theoretically_ work this out, but due to
+	 * the nature of floating point arithmetic and rounding errors, we need to
+	 * crunch the numbers the hard way */
+	// size_t new_width = total_width * pix_per_deg;
+	// size_t new_height = total_height * pix_per_deg;
+
+	size_t new_height = 0;
+	size_t new_width = 0;
+	for ( size_t i = 0; i < fc; i++ ) {
+		double north_offset = max_north - tiles[i].max_north;
+		double west_offset = max_west - tiles[i].max_west >= 0 ? max_west - tiles[i].max_west : max_west + (360 - tiles[i].max_west);
+		size_t north_pixel_offset = north_offset * tiles[i].ppdy;
+		size_t west_pixel_offset = west_offset * tiles[i].ppdx;
+
+		if ( west_pixel_offset + tiles[i].width > new_width )
+			new_width = west_pixel_offset + tiles[i].width;
+		if ( north_pixel_offset + tiles[i].height > new_height )
+			new_height = north_pixel_offset + tiles[i].height;
+	}
+	size_t new_tile_alloc = new_width * new_height;
+
+	short * new_tile = (short*) calloc( new_tile_alloc, sizeof(short) );
+	if ( new_tile == NULL ){
+		free(tiles);
+		return ENOMEM;
+	}
+	if (debug)
+		fprintf(stderr,"Lidar tile dimensions w:%lf(%zu) h:%lf(%zu)\n", total_width, new_width, total_height, new_height);
+
+	/* ...If we wanted a value other than sea level here, we would need to initialize the array... */
+
+	/* Fill out the array one tile at a time */
+	for (size_t i = 0; i< fc; i++) {
+		double north_offset = max_north - tiles[i].max_north;
+		double west_offset = max_west - tiles[i].max_west >= 0 ? max_west - tiles[i].max_west : max_west + (360 - tiles[i].max_west);
+		size_t north_pixel_offset = north_offset * tiles[i].ppdy;
+		size_t west_pixel_offset = west_offset * tiles[i].ppdx;
 
 
-		fgets(line, 255, fd); // NODATA
-		pos = ftell(fd);
+		if (debug) {
+			fprintf(stderr,"mn: %lf mw:%lf globals: %lf %lf\n", tiles[i].max_north, tiles[i].max_west, max_north, max_west);
+			fprintf(stderr,"Offset n:%zu(%lf) w:%zu(%lf)\n", north_pixel_offset, north_offset, west_pixel_offset, west_offset);
+			fprintf(stderr,"Height: %d\n", tiles[i].height);
+		}
 
-		// tile 0 [x| ]
-		if (debug)
-			fprintf(stderr, "readLIDAR(fd,%d,%d,%d,%.4f,%.4f,%.4f,%.4f)\n", height, width, indx, yur, xur, yll, xll);
+		/* Copy it row-by-row from the tile */
+		for (size_t h = 0; h < tiles[i].height; h++) {
+			register short *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset];
+			register short *src_addr = &tiles[i].data[h*tiles[i].width];
+			// Check if we might overflow
+			if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc || dest_addr < new_tile ){
+				if (debug)
+					fprintf(stderr, "Overflow %zu\n",i);
+				continue;
+			}
+			memcpy( dest_addr, src_addr, tiles[i].width * sizeof(short) );
+		}
+	}
 
-		readLIDAR(fd, height, width, indx, yur, xur, yll, xll);
+	MAXPAGES = 1;
+	IPPD = MAX(new_width,new_height);
+	if(debug){
+		fprintf(stderr,"Setting IPPD to %d\n",IPPD);
+		fflush(stderr);
+	}
+	ARRAYSIZE = (MAXPAGES * IPPD) + 50;
+	do_allocs();
+
+	/* Load the data into the global dem array */
+	dem[0].max_north = max_north;
+	dem[0].min_west = min_west;
+	dem[0].min_north = min_north;
+	dem[0].max_west = max_west;
+	dem[0].max_el = max_elevation;
+	dem[0].min_el = min_elevation;
+
+	/* 
+	 * Copy the lidar tile data into the dem array. The dem array is rotated
+	 * 90 degrees (christ knows why...)
+	 */
+	int y = new_height-1;
+	for (size_t h = 0; h < new_height; h++, y--) {
+		int x = new_width-1;
+		for (size_t w = 0; w < new_width; w++, x--) {
+			dem[0].data[y][x] = new_tile[h*new_width + w];
+			dem[0].signal[y][x] = 0;
+			dem[0].mask[y][x] = 0;
+		}
+	}
 
-		fclose(fd);
-		if (debug)
-			fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height);
-		indx++;
-	} // filename(s)
-	IPPD=width;
 	ippd=IPPD;
-	height = (unsigned)((max_north-min_north) / cellsize);
-	width = (unsigned)((max_west-min_west) / cellsize);
+	height = new_height;
+	width = new_width;
+
+	if (debug)
+		fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height);
 
 	if (debug)
 		fprintf(stderr, "fc %d WIDTH %d HEIGHT %d ippd %d minN %.5f maxN %.5f minW %.5f maxW %.5f avgCellsize %.5f\n", fc, width, height, ippd,min_north,max_north,min_west,max_west,avgCellsize);
+
+cleanup:
+
+	if ( tiles != NULL ) {
+		for (size_t i = 0; i < fc; i++) {
+			tile_destroy(&tiles[i]);
+		}
+	}
+	free(tiles);
+
 	return 0;
 }
 
diff --git a/inputs.hh b/inputs.hh
index c946286..6d7d353 100644
--- a/inputs.hh
+++ b/inputs.hh
@@ -3,6 +3,10 @@
 
 #include "common.h"
 
+/* Resample input tiles to new resolution */
+int resample_data(int scaling_factor);
+int resize_data(int resolution);
+
 int LoadSDF_SDF(char *name, int winfiles);
 int LoadSDF(char *name, int winfiles);
 int LoadPAT(char *az_filename, char *el_filename);
@@ -11,7 +15,7 @@ int LoadLossColors(struct site xmtr);
 int LoadDBMColors(struct site xmtr);
 int LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat);
 int LoadUDT(char *filename);
-int loadLIDAR(char *filename);
+int loadLIDAR(char *filename, int resample);
 int loadClutter(char *filename, double radius, struct site tx);
 
 static const char AZ_FILE_SUFFIX[] = ".az";
diff --git a/main.cc b/main.cc
index 7634848..cc19762 100644
--- a/main.cc
+++ b/main.cc
@@ -53,11 +53,11 @@ double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd,
 int ippd, mpi,
     max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold,
     pred, pblue, pgreen, ter, multiplier = 256, debug = 0, loops = 100, jgets =
-    0, MAXRAD, hottest = 10, height, width;
+    0, MAXRAD, hottest = 10, height, width, resample = 0;
 
 unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0;
 
-bool to_stdout = false;
+bool to_stdout = false, cropping = false;
 
 __thread double *elev;
 __thread struct path path;
@@ -1097,6 +1097,7 @@ int main(int argc, char *argv[])
 		fprintf(stdout,	"     -tercon Terrain conductivity 0.01-0.0001 (optional)\n");
 		fprintf(stdout, "     -cl Climate code 1-6 (optional)\n");
 		fprintf(stdout, "     -rel Reliability for ITM model 50 to 99 (optional)\n");
+		fprintf(stdout, "     -resample Resample Lidar input to specified resolution in meters (optional)\n");
 		fprintf(stdout, "Output:\n");
 		fprintf(stdout,	"     -dbm Plot Rxd signal power instead of field strength\n");
 		fprintf(stdout, "     -rt Rx Threshold (dB / dBm / dBuV/m)\n");
@@ -1145,6 +1146,7 @@ int main(int argc, char *argv[])
 	max_txsites = 30;
 	fzone_clearance = 0.6;
 	contour_threshold = 0;
+	resample = 0;
 
 	ano_filename[0] = 0;
 	earthradius = EARTHRADIUS;
@@ -1322,6 +1324,17 @@ int main(int argc, char *argv[])
 			}
 		}
 
+		if (strcmp(argv[x], "-resample") == 0) {
+			z = x + 1;
+
+			if(!lidar){
+				fprintf(stderr, "Error, this should only be used with LIDAR tiles.\n");
+				return -1;
+			}
+
+			sscanf(argv[z], "%d", &resample);
+		}
+
 		if (strcmp(argv[x], "-lat") == 0) {
 			z = x + 1;
 
@@ -1687,10 +1700,11 @@ int main(int argc, char *argv[])
 
 	/* Load the required tiles */
 	if(lidar){
-		if( (result = loadLIDAR(lidar_tiles)) != 0 ){
+		if( (result = loadLIDAR(lidar_tiles, resample)) != 0 ){
 			fprintf(stderr, "Couldn't find one or more of the "
 				"lidar files. Please ensure their paths are "
 				"correct and try again.\n");
+			fprintf(stderr, "Error %d: %s\n", result, strerror(result));
 			exit(result);
 		}
 
@@ -1710,9 +1724,9 @@ int main(int argc, char *argv[])
 
 	}else{
 		// DEM first
-                if(debug){
-                        fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
-                }
+		if(debug){
+			fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
+		}
 
 		max_lon-=3;
 
@@ -1867,6 +1881,7 @@ int main(int argc, char *argv[])
 				}
 			}
 
+#ifdef CROPPING
 			if(max_range<=100){
 				// CROPPING. croplat assigned in propPathLoss()
 				max_north=cropLat; // MAX(path.lat[y])
@@ -1881,12 +1896,14 @@ int main(int argc, char *argv[])
 				cropLon-=tx_site[0].lon;
 				width=(int)((cropLon*ppd)*2);
 				height=(int)((cropLat*ppd)*2);
+				cropping = true;
 
 				if(width>3600*10){
 					fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon);
 					return 0;
 				}
 			}
+#endif
 
 			// Write bitmap
 			if (LR.erp == 0.0)
@@ -1908,21 +1925,19 @@ int main(int argc, char *argv[])
 					tx_site[0].lon *= -1;
 		}
 		if (tx_site[0].lon < -180.0){
-					tx_site[0].lon += 360;
+			tx_site[0].lon += 360;
 		}
 
-		if (propmodel == 2 || max_range > 100) {
-			// No croppping because this is LOS
-			fprintf(stderr, "|%.6f", max_north);
-			fprintf(stderr, "|%.6f", east);
-			fprintf(stderr, "|%.6f", min_north);
-			fprintf(stderr, "|%.6f|",west);
-		}else{
-			// Cropped EPSG4326 coordinates
+		if (cropping) {
 			fprintf(stderr, "|%.6f", tx_site[0].lat+cropLat);
 			fprintf(stderr, "|%.6f", tx_site[0].lon+cropLon);
 			fprintf(stderr, "|%.6f", tx_site[0].lat-cropLat);
 			fprintf(stderr, "|%.6f|",tx_site[0].lon-cropLon);
+		}else{
+			fprintf(stderr, "|%.6f", max_north);
+			fprintf(stderr, "|%.6f", east);
+			fprintf(stderr, "|%.6f", min_north);
+			fprintf(stderr, "|%.6f|",west);
 		}
 		fprintf(stderr, "\n");
 
diff --git a/outputs.cc b/outputs.cc
index 5d79961..d132849 100644
--- a/outputs.cc
+++ b/outputs.cc
@@ -35,7 +35,7 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml,
 	int success;
 
 	if( (success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0 ){
-		fprintf(stderr,"Error initializing image\n");
+		fprintf(stderr,"Error initializing image: %s\n", strerror(success));
 		exit(success);
 	}
 
@@ -278,7 +278,7 @@ int DoSigStr(char *filename, unsigned char geo, unsigned char kml,
 	int success;
 
 	if((success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0){
-		fprintf(stderr,"Error initializing image\n");
+		fprintf(stderr,"Error initializing image: %s\n", strerror(success));
 		exit(success);
 	}
 
@@ -529,7 +529,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml,
 	int success;
 
 	if( (success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0 ){
-		fprintf(stderr,"Error initializing image\n");
+		fprintf(stderr,"Error initializing image: %s\n", strerror(success));
 		exit(success);
 	}
 
@@ -775,7 +775,7 @@ void DoLOS(char *filename, unsigned char geo, unsigned char kml,
 	int success;
 
 	if((success = image_init(&ctx, width, (kml ? height : height + 30), IMAGE_RGB, IMAGE_DEFAULT)) != 0){
-		fprintf(stderr,"Error initializing image\n");
+		fprintf(stderr,"Error initializing image: %s\n", strerror(success));
 		exit(success);
 	}
 
diff --git a/tiles.cc b/tiles.cc
new file mode 100644
index 0000000..f81bcb7
--- /dev/null
+++ b/tiles.cc
@@ -0,0 +1,245 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <math.h>
+#include "tiles.hh"
+#include "common.h"
+
+#define MAX_LINE 25000
+
+/* Computes the distance between two long/lat points */
+double haversine_formula(double th1, double ph1, double th2, double ph2)
+{
+	#define TO_RAD (3.1415926536 / 180)
+	int R = 6371;
+	double dx, dy, dz;
+	ph1 -= ph2;
+	ph1 *= TO_RAD, th1 *= TO_RAD, th2 *= TO_RAD;
+ 
+	dz = sin(th1) - sin(th2);
+	dx = cos(ph1) * cos(th1) - cos(th2);
+	dy = sin(ph1) * cos(th1);
+	return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * R;
+}
+
+int tile_load_lidar(tile_t *tile, char *filename){
+	FILE *fd;
+	char line[MAX_LINE];
+	short nextval;
+	char *pch;
+
+	/* Clear the tile data */
+	memset(tile, 0x00, sizeof(tile_t));
+
+	/* Open the file handle and return on error */
+	if ( (fd = fopen(filename,"r")) == NULL )
+		return errno;
+
+	/* This is where we read the header data */
+	/* The string is split for readability but is parsed as a block */
+	if( fscanf(fd,"%*s %d\n" "%*s %d\n" "%*s %lf\n" "%*s %lf\n" "%*s %lf\n" "%*s %d\n",&tile->width,&tile->height,&tile->xll,&tile->yll,&tile->cellsize,&tile->nodata) != 6 ){
+		fclose(fd);
+		return -1;
+	}
+
+	tile->datastart = ftell(fd);
+
+	if(debug){
+		fprintf(stderr,"w:%d h:%d s:%lf\n", tile->width, tile->height, tile->cellsize);
+		fflush(stderr);
+	}
+
+	/* Set the filename */
+	tile->filename = strdup(filename);
+
+	/* Perform xur calcs */
+	// Degrees with GDAL option: -co "FORCE_CELLSIZE=YES"
+	tile->xur = tile->xll+(tile->cellsize*tile->width);
+	tile->yur = tile->yll+(tile->cellsize*tile->height);
+
+	if (tile->xur > eastoffset)
+		eastoffset = tile->xur;
+	if (tile->xll < westoffset)
+		westoffset = tile->xll;
+
+	// if (debug)
+	// 	fprintf(stderr,"%d, %d, %.7f, %.7f, %.7f, %.7f, %.7f\n",width,height,xll,yll,cellsize,yur,xur);
+
+	// Greenwich straddling hack
+	// if (tile->xll <= 0 && tile->xur > 0) {
+	// 	tile->xll = (tile->xur - tile->xll); // full width
+	// 	tile->xur = 0.0; // budge it along so it's west of greenwich
+	// 	delta = eastoffset; // add to Tx longitude later
+	// } else {
+		// Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;)
+		if (tile->xll >= 0)
+			tile->xll = 360-tile->xll;
+		if(tile->xur >= 0)
+			tile->xur = 360-tile->xur;
+		if(tile->xll < 0)
+			tile->xll = tile->xll * -1;
+		if(tile->xur < 0)
+			tile->xur = tile->xur * -1;
+	// }
+
+	if (debug)
+		fprintf(stderr, "POST yll %.7f yur %.7f xur %.7f xll %.7f delta %.6f\n", tile->yll, tile->yur, tile->xur, tile->xll, delta);
+
+	/* Read the actual tile data */
+	/* Allocate the array for the lidar data */
+	if ( (tile->data = (short*) calloc(tile->width * tile->height, sizeof(short))) == NULL ) {
+		fclose(fd);
+		free(tile->filename);
+		return ENOMEM;
+	}
+
+	size_t loaded = 0;
+	for (size_t h = 0; h < tile->height; h++) {
+		if (fgets(line, MAX_LINE, fd) != NULL) {
+			pch = strtok(line, " "); // split line into values
+			for (size_t w = 0; w < tile->width && pch != NULL; w++) {
+				/* If the data is less than a *magic* minimum, normalize it to zero */
+				nextval = atoi(pch);
+				if (nextval <= -9999)
+					nextval = 0;
+				tile->data[h*tile->width + w] = nextval;
+				loaded++;
+				if ( nextval > tile->max_el )
+					tile->max_el = nextval;
+				if ( nextval < tile->min_el )
+					tile->min_el = nextval;
+				pch = strtok(NULL, " ");
+			}//while
+		} else {
+			fprintf(stderr, "LIDAR error @ h %zu file %s\n", h, filename);
+		}//if
+	}
+
+	double current_res_km = haversine_formula(tile->max_north, tile->max_west, tile->max_north, tile->min_west);
+	tile->resolution = (int) ceil((current_res_km/MAX(tile->width,tile->height))*1000);
+
+	tile->width_deg = tile->max_west - tile->min_west >= 0 ? tile->max_west - tile->min_west : tile->max_west + (360 - tile->min_west);
+	tile->height_deg = tile->max_north - tile->min_north;
+
+	tile->ppdx = tile->width / tile->width_deg;
+	tile->ppdy = tile->height / tile->height_deg;
+
+	if (debug)
+		fprintf(stderr,"Pixels loaded: %zu/%d (PPD %dx%d)\n", loaded, tile->width*tile->height, tile->ppdx, tile->ppdy);
+
+	/* All done, close the LIDAR file */
+	fclose(fd);
+
+	return 0;
+}
+
+/*
+ * tile_rescale
+ * This is used to resample tile data. It is particularly designed for
+ * use with LIDAR tiles where the resolution can be anything up to 2m.
+ * This function is capable of merging neighbouring pixel values
+ * The scaling factor is the distance to merge pixels.
+ * NOTE: This means that new resolutions can only increment in multiples of the original
+ * (ie 2m LIDAR can be 4/6/8/... and 20m can be 40/60)
+ */
+int tile_rescale(tile_t *tile, float scale){
+	short *new_data;
+	size_t skip_count = 1;
+	size_t copy_count = 1;
+
+	if (scale == 1) {
+		return 0;	
+	}
+
+	size_t new_height = tile->height * scale;
+	size_t new_width = tile->width * scale;
+
+	/* Allocate the array for the lidar data */
+	if ( (new_data = (short*) calloc(new_height * new_width, sizeof(short))) == NULL ) {
+		return ENOMEM;
+	}
+
+	tile->max_el = -32768;
+	tile->min_el = 32768;
+
+	/* Making the tile data smaller */
+	if (scale < 1) {
+		skip_count = 1 / scale;
+	} else {
+		copy_count = (size_t) scale;
+	}
+
+	if (debug)
+		fprintf(stderr,"Resampling tile:\n\tOld %zux%zu. New %zux%zu\n\tScale %f Skip %zu Copy %zu\n", tile->width, tile->height, new_width, new_height, scale, skip_count, copy_count);
+
+	/* Nearest neighbour normalization. For each subsample of the original, simply
+	 * assign the value in the top left to the new pixel 
+	 * SOURCE: X / Y
+	 * DEST:   I / J */
+
+	for (size_t y = 0, j = 0;
+			y < tile->height && j < new_height;
+			y += skip_count, j += copy_count){
+
+		for (size_t x = 0, i = 0;
+				x < tile->width && i < new_width;
+				x += skip_count, i += copy_count) {
+		
+			/* These are for scaling up the data */
+			for (size_t copy_y = 0; copy_y < copy_count; copy_y++) {
+				for (size_t copy_x = 0; copy_x < copy_count; copy_x++) {
+					size_t new_j = j + copy_y;
+					size_t new_i = i + copy_x;
+					/* Do the copy */
+					new_data[ new_j * new_width + new_i ] = tile->data[y * tile->width + x];
+				}
+			}
+			/* Update local min / max values */
+			if (tile->data[y * tile->width + x] > tile->max_el)
+				tile->max_el = tile->data[y * tile->width + x];
+			if (tile->data[y * tile->width + x] < tile->min_el)
+				tile->min_el = tile->data[y * tile->width + x];
+		}
+	}
+
+	/* Update the date in the tile */
+	free(tile->data);
+	tile->data = new_data;
+
+	/* Update the height and width values */
+	tile->height = new_height;
+	tile->width = new_width;
+	tile->resolution *= scale;
+	tile->ppdy *= scale;
+	tile->ppdx *= scale;
+	tile->width_deg *= scale;
+	tile->height_deg *= scale;
+
+	return 0;
+}
+
+/*
+ * tile_resize
+ * This function works in conjuntion with resample_data. It takes a
+ * resolution value in meters as its argument. It then calculates the
+ * nearest (via averaging) resample value and calls resample_data
+ */
+int tile_resize(tile_t* tile, int resolution){
+	double current_res_km = haversine_formula(tile->max_north, tile->max_west, tile->max_north, tile->min_west);
+	int current_res = (int) ceil((current_res_km/IPPD)*1000);
+	float scaling_factor = resolution / current_res;
+	if (debug)
+		fprintf(stderr, "Resampling: Current %dm Desired %dm Scale %d\n", current_res, resolution, scaling_factor);
+	return tile_rescale(tile, scaling_factor);
+}
+
+/*
+ * tile_destroy
+ * This function simply destroys any data associated with a tile
+ */
+void tile_destroy(tile_t* tile){
+	if (tile->data != NULL)
+		free(tile->data);
+}
+
diff --git a/tiles.hh b/tiles.hh
new file mode 100644
index 0000000..a2d25f4
--- /dev/null
+++ b/tiles.hh
@@ -0,0 +1,47 @@
+#ifndef _TILES_HH_
+#define _TILES_HH_
+
+typedef struct _tile_t{
+	char	*filename;
+	union{
+		int	cols;
+		int	width;
+	};
+	union{
+		int	rows;
+		int	height;
+	};
+	union{
+		double	xll;
+		double	max_west;
+	};
+	union{
+		double	yll;
+		double	min_north;
+	};
+	union{
+		double	xur;
+		double	min_west;
+	};
+	union{
+		double	yur;
+		double	max_north;
+	};
+	double	cellsize;
+	long long datastart;
+	short	nodata;
+	short 	max_el;
+	short	min_el;
+	short	*data;
+	int		resolution;
+	double	width_deg;
+	double	height_deg;
+	int		ppdx;
+	int		ppdy;
+} tile_t, *ptile_t;
+
+int tile_load_lidar(tile_t*, char *);
+int tile_rescale(tile_t *, float);
+void tile_destroy(tile_t *);
+
+#endif
\ No newline at end of file