Merge pull request #6 from emilyselwood/test-branch

updates to make stats functions work

Author: emilyselwood

datacube_utilities/dc_mosaic.py

@@ -61,7 +61,7 @@ def convert_to_dtype(data, dtype):
 Compositing Functions
 """
 
-def create_mosaic(dataset_in, clean_mask=None, no_data=-9999, dtype=None, intermediate_product=None, **kwargs):
+def create_mosaic(dataset_in, clean_mask=None, no_data=float('nan'), dtype=None, intermediate_product=None, **kwargs):
     """
     Creates a most-recent-to-oldest mosaic of the input dataset.
 
@@ -125,7 +125,7 @@ def create_mosaic(dataset_in, clean_mask=None, no_data=-9999, dtype=None, interm
     dataset_out = restore_or_convert_dtypes(dtype, band_list, dataset_in_dtypes, dataset_out, no_data)
     return dataset_out
 
-def create_mean_mosaic(dataset_in, clean_mask=None, no_data=-9999, dtype=None, **kwargs):
+def create_mean_mosaic(dataset_in, clean_mask=None, no_data=float('nan'), dtype=None, **kwargs):
     """
     Method for calculating the mean pixel value for a given dataset.
 
@@ -215,7 +215,7 @@ def create_median_mosaic(dataset_in, clean_mask=None, no_data=float('nan'), dtyp
     return dataset_out
 
 
-def create_max_ndvi_mosaic(dataset_in, clean_mask=None, no_data=-9999, dtype=None, intermediate_product=None, **kwargs):
+def create_max_ndvi_mosaic(dataset_in, clean_mask=None, no_data=float('nan'), dtype=None, intermediate_product=None, **kwargs):
     """
     Method for calculating the pixel value for the max ndvi value.
 
@@ -279,7 +279,7 @@ def create_max_ndvi_mosaic(dataset_in, clean_mask=None, no_data=-9999, dtype=Non
     return dataset_out
 
 
-def create_min_ndvi_mosaic(dataset_in, clean_mask=None, no_data=[-9999, 0, np.float('nan')], dtype=None, intermediate_product=None, **kwargs):
+def create_min_ndvi_mosaic(dataset_in, clean_mask=None, no_data=float('nan'), dtype=None, intermediate_product=None, **kwargs):
     """
     Method for calculating the pixel value for the min ndvi value.
 
@@ -525,12 +525,12 @@ def restore_or_convert_dtypes(dtype_for_all, band_list, dataset_in_dtypes, datas
     if dtype_for_all is not None:
         # Integer types can't represent nan.
         if np.issubdtype(dtype_for_all, np.integer): # This also works for Python int type.
-            utilities.nan_to_num(dataset_out, no_data)
+            dataset_out.fillna(no_data)
         convert_to_dtype(dataset_out, dtype_for_all)
     else:  # Restore dtypes to state before masking.
         for band in band_list:
             band_dtype = dataset_in_dtypes[band]
             if np.issubdtype(band_dtype, np.integer):
-                utilities.nan_to_num(dataset_out[band], no_data)
+                dataset_out[band].fillna(no_data)
             dataset_out[band] = dataset_out[band].astype(band_dtype)
     return dataset_out

datacube_utilities/dc_utilities.py

@@ -182,11 +182,7 @@ def nan_to_num(data, number):
     ----------
     data: xarray.Dataset or xarray.DataArray
     """
-    if isinstance(data, xr.Dataset):
-        for key in list(data.data_vars):
-            data[key].values[np.isnan(data[key].values)] = number
-    elif isinstance(data, xr.DataArray):
-        data.values[np.isnan(data.values)] = number
+    data.fillna(number)
 
 
 def clear_attrs(dataset):

datacube_utilities/query.py

@@ -4,16 +4,32 @@
 from datacube_utilities.createAOI import create_lat_lon
 
 
-def create_base_query(aoi, res, output_projection, aoi_crs, dask_chunks):
+def create_base_query(aoi, res, output_crs, aoi_crs, dask_chunks, cube_crs="EPSG:3460"):
+    """
+    create_base_query sets up the basic data cube query and makes sure the AOI is in the right CRS
+
+    Parameters
+    ----------
+    aoi: WKT formatted string.
+    res: Int Resolution of the output images
+    output_crs: Resulting projection string EPSG code.
+    aoi_crs: Projection string EPSG code of the AOI . Does not have to be the same as the output.
+    dask_chunks: Number of dask chunks to split the resulting data into.
+    cube_crs: The Projection string EPSG code used by the data cube.
+
+    Returns
+    -------
+    query: A query object ready to be fed into dc.load
+    """
     lat_extents, lon_extents = create_lat_lon(aoi)
-    inProj = Proj("+init=EPSG:4326")
-    outProj = Proj("+init=EPSG:3460")
+    in_proj = Proj("+init=" + aoi_crs)
+    out_proj = Proj("+init=" + cube_crs)
 
     min_lat, max_lat = lat_extents
     min_lon, max_lon = lon_extents
-
-    x_A, y_A = transform(inProj, outProj, min_lon, min_lat)
-    x_B, y_B = transform(inProj, outProj, max_lon, max_lat)
+    
+    x_A, y_A = transform(in_proj, out_proj, min_lon, min_lat)
+    x_B, y_B = transform(in_proj, out_proj, max_lon, max_lat)
 
     lat_range = (y_A, y_B)
     lon_range = (x_A, x_B)
@@ -23,10 +39,10 @@ def create_base_query(aoi, res, output_projection, aoi_crs, dask_chunks):
     query = {
         "y": lat_range,
         "x": lon_range,
-        "output_crs": output_projection,
+        "output_crs": output_crs,
         "resolution": resolution,
         "dask_chunks": dask_chunks,
-        "crs": aoi_crs,
+        "crs": cube_crs,
     }
     return query