diff --git a/README.md b/README.md
index a6100f3..1aca843 100644
--- a/README.md
+++ b/README.md
@@ -140,7 +140,7 @@ Example Usage: Other functions
 	NDV, xsize, ysize, GeoT, Projection, DataType = gr.get_geo_info(raster)
 	
 	# Load raster
-	data = load_tiff(raster)
+	data = gr.load_tiff(raster)
 	
 	# Find location of point (x,y) on raster, e.g. to extract info at that location
 	col, row = gr.map_pixel(x,y,GeoT[1],GeoT[-1], GeoT[0],GeoT[3])
@@ -150,7 +150,7 @@ Example Usage: Other functions
 	gr.aggregate(data,NDV,(10,10))
 	
 	# Align two rasters
-	data2 = load_tiff(raster2)
+	data2 = gr.load_tiff(raster2)
 	(alignedraster_o, alignedraster_a, GeoT_a) = gr.align_rasters(raster, raster2, how=np.mean)
 	
 	# Create GeoRaster
@@ -158,7 +158,7 @@ Example Usage: Other functions
 	
 	# Load another raster
 	NDV, xsize, ysize, GeoT, Projection, DataType = gr.get_geo_info(raster2)
-	data = load_tiff(raster2)
+	data = gr.load_tiff(raster2)
 	B=gr.GeoRaster(data2, GeoT, nodata_value=NDV)
 	
 	# Plot Raster
diff --git a/README.rst b/README.rst
index e648d3d..db0f013 100644
--- a/README.rst
+++ b/README.rst
@@ -173,7 +173,7 @@ Example Usage: Other functions
     NDV, xsize, ysize, GeoT, Projection, DataType = gr.get_geo_info(raster)
     
     # Load raster
-    data = load_tiff(raster)
+    data = gr.load_tiff(raster)
        
     # Find location of point (x,y) on raster, e.g. to extract info at that location
     col, row = gr.map_pixel(x,y,GeoT[1],GeoT[-1], GeoT[0],GeoT[3])
@@ -183,7 +183,7 @@ Example Usage: Other functions
     gr.aggregate(data,NDV,(10,10))
     
     # Align two rasters
-    data2 = load_tiff(raster2)
+    data2 = gr.load_tiff(raster2)
     (alignedraster_o, alignedraster_a, GeoT_a) = gr.align_rasters(raster, raster2, how=np.mean)
     
     # Create GeoRaster
@@ -191,7 +191,7 @@ Example Usage: Other functions
 
     # Load another raster
     NDV, xsize, ysize, GeoT, Projection, DataType = gr.get_geo_info(raster2)
-    data = load_tiff(raster2)
+    data = gr.load_tiff(raster2)
     B=gr.GeoRaster(data2, GeoT, nodata_value=NDV)
     
     # Plot Raster
diff --git a/georasters/georasters.py b/georasters/georasters.py
index 2cf9d7d..6405f99 100755
--- a/georasters/georasters.py
+++ b/georasters/georasters.py
@@ -114,24 +114,25 @@ def aggregate(raster, ndv, block_size):
     return raster2
 
 # Function to write a new file.
-def create_geotiff(name, Array, driver, ndv, xsize, ysize, geot, projection, datatype, band=1):
+def create_geotiff(name, Array, driver, ndv, xsize, ysize, geot, projection, datatype, band=1, **kwargs):
     '''
-    Creates new geotiff from array
+    Creates new geotiff from array.
     '''
     if isinstance(datatype, np.int) == False:
         if datatype.startswith('gdal.GDT_') == False:
             datatype = eval('gdal.GDT_'+datatype)
     newfilename = name+'.tif'
     # Set nans to the original No Data Value
-    Array[np.isnan(Array)] = ndv
+    if len(Array[np.isnan(Array)]) > 0:
+        Array[np.isnan(Array)] = ndv
     # Set up the dataset
-    DataSet = driver.Create(newfilename, xsize, ysize, 1, datatype)
+    DataSet = driver.Create(newfilename, xsize, ysize, 1, datatype, **kwargs)
     # the '1' is for band 1.
     DataSet.SetGeoTransform(geot)
     DataSet.SetProjection(projection.ExportToWkt())
     # Write the array
     DataSet.GetRasterBand(band).WriteArray(Array)
-    DataSet.GetRasterBand(band).SetNoDataValue(ndv)
+    DataSet.GetRasterBand(band).SetNoDataValue(np.float64(ndv))
     return newfilename
 
 # Function to aggregate and align rasters
@@ -455,7 +456,7 @@ def copy(self):
         return GeoRaster(self.raster.copy(), self.geot, nodata_value=self.nodata_value,
                          projection=self.projection, datatype=self.datatype)
 
-    def to_tiff(self, filename):
+    def to_tiff(self, filename, **kwargs):
         '''
         geo.to_tiff(filename)
 
@@ -464,6 +465,10 @@ def to_tiff(self, filename):
         If GeoRaster does not have datatype, then it tries to assign a type.
         You can assign the type yourself by setting
          geo.datatype = 'gdal.GDT_'+type
+
+        kwargs are passed to GeoTiff driver, and may include GDAL options like compression:
+
+            ds.to_tiff(filename, options=['COMPRESS=LZMA'])
         '''
         if self.datatype is None:
             self.datatype = gdal_array.NumericTypeCodeToGDALTypeCode(self.raster.data.dtype)
@@ -474,9 +479,11 @@ def to_tiff(self, filename):
                 else:
                     self.raster = self.raster.astype(np.float64)
                     self.datatype = gdal_array.NumericTypeCodeToGDALTypeCode(self.raster.data.dtype)
-        self.raster.data[self.raster.mask] = self.nodata_value
+        
+        if len(self.raster.data[self.raster.mask]) > 0:
+            self.raster.data[self.raster.mask] = self.nodata_value
         create_geotiff(filename, self.raster, gdal.GetDriverByName('GTiff'), self.nodata_value,
-                       self.shape[1], self.shape[0], self.geot, self.projection, self.datatype)
+                       self.shape[1], self.shape[0], self.geot, self.projection, self.datatype, **kwargs)
 
     def to_pandas(self):
         """