Repository: tnigon/hs_process (language: python)

Issue:
Crop image issue
Hi, 
I am cropping UAV hyperspectral imagery with plot boundary files, there is one issue with this line of code:
```
array_crop, metadata = my_spatial_mod.crop_single(plot_id=1)


KeyError                                  Traceback (most recent call last)
<ipython-input-43-458553a456fd> in <module>
----> 1 array_crop, metadata = my_spatial_mod.crop_single(plot_id=1)

~\Anaconda3\envs\spec\lib\site-packages\hs_process\spatial_mod.py in crop_single(self, pix_e_ul, pix_n_ul, crop_e_pix, crop_n_pix, crop_e_m, crop_n_m, buf_e_pix, buf_n_pix, buf_e_m, buf_n_m, spyfile, plot_id, gdf, name_append)
    909                     "".format(pix_e_ul, pix_n_ul, pix_e_lr, pix_n_lr))
    910         # If "..crop_single" is already included in the history, remove it
--> 911         idx_remove = metadata['history'].find(
    912                 ' -> hs_process.crop_single[<')
    913         if idx_remove == -1:

KeyError: 'history'
```

Any idea how to solve the issue? 
Thanks. 


Discussion / hints:
The error tells me that the image does not have a "history" key in its metadata. A workaround might be to simply add a "history" tag to the metadata, even if it is blank or just has a dummy variable. Can you provide more of your code? What is the image type you have loaded into ``spatial_mod``? Can you share the contents of your metadata? I'll try to help you out if I can. Thanks!
> The error tells me that the image does not have a "history" key in its metadata. A workaround might be to simply add a "history" tag to the metadata, even if it is blank or just has a dummy variable. Can you provide more of your code? What is the image type you have loaded into `spatial_mod`? Can you share the contents of your metadata? I'll try to help you out if I can. Thanks!

My code is about the same as section 6. Tutorial: spatial_mod:
The metadata seems not defined. But I am not quite sure how to create the metadata. 
```
import geopandas as gpd
import os
from hs_process import hsio
from hs_process import spatial_mod

fname_in = r'E:\hyperspectral_imagery1.hdr'
fname_gdf = r'E:\plotboundary.geojson'
gdf = gpd.read_file(fname_gdf)
io = hsio(fname_in)
my_spatial_mod = spatial_mod(io.spyfile)
dir_out = os.path.join(io.base_dir, 'spatial_mod', 'crop_many_gdf')
if not os.path.isdir(os.path.join(io.base_dir, 'spatial_mod')):  # create a new folder named "spatial_mod" if it does not exist
    os.mkdir(os.path.join(io.base_dir, 'spatial_mod'))
if not os.path.isdir(dir_out):  # create a new folder named "crop_many_gdf" if it does not exist
    os.mkdir(dir_out)

name_append = '-crop-many-gdf'
df_plots = my_spatial_mod.crop_many_gdf(spyfile=io.spyfile, gdf=gdf)
array_crop, metadata = my_spatial_mod.crop_single(plot_id=1)
```
Thanks for the information. Can you open up the .hdr file (``E:\hyperspectral_imagery1.hdr``) with a text editor and look at the contents? Here is a sample of the image I used in the example/tutorial:

```
ENVI
samples = 3
lines = 3
bands = 240
header offset = 0
file type = ENVI Standard
data type = 4
interleave = bip
byte order = 0
reflectance scale factor = 1.0
map info = {UTM, 1.0, 1.0, 441371.00707299996, 4855941.7317699995, 0.04, 0.04, 15, T, WGS-84, units  meters, rotation  0.000}
history = input cube -> Radiance Conversion<SpecGroup>[<SpecFilename label:'Imager Calibration' value:J:\\\\Shared Work\\\\Data\\\\PikaImagery2_Radiance\\\\Pika II 100106_11 Radiometric Calibration September 2016\\\\8mm\\\\PikaII_100106_11_8mm_RadiometricCalibration_Sept2016.icp>, <SpecCube label:'Dark Noise Cube' value:None>, <SpecBool label:'Auto Remove Dark Noise?' value:True>, <SpecBool label:'Return floating point?' value:False>] -> Georectify Airborne Datacube<SpecGroup>[[<SpecBool label:'Use Flat Earth Altitude?' value:False>, <SpecFloatText label:'Flat Earth Alt (m)' value:340>, <SpecFloatText label:'Field of View (deg)' value:33.0>, <SpecInfo label:'Calculated Resolution (m)' value:Not Available in Batch Mode>, <SpecFloatText label:'Map Resolution (m)' value:0.04>, <SpecBool label:'Select DEM File' value:False>], [<SpecFloatText label:'Sync Offset (s)' value:-0.05>, <SpecFloatText label:'Imager Roll Offset (deg)' value:0>, <SpecFloatText label:'Imager Pitch Offset (deg)' value:2.25>, <SpecFloatText label:'Imager Heading Offset (deg)' value:0>], [<SpecBool label:'Crop Cube to Start/End Lines?' value:False>, <SpecIntText label:'Start Line' value:0>, <SpecIntText label:'End Line' value:1906>], [<SpecBool label:'Select Folder to Save Output' value:False>, <SpecBool label:'Save Products in Source Folder?' value:False>, <SpecBool label:'Generate Full Datacube?' value:True>, <SpecBool label:'Generate KML of Render?' value:True>, <SpecBool label:'Generate GeoTIFF Render?' value:True>, <SpecBool label:'Generate Swath Outline?' value:True>, <SpecChoice label:'Interpolate Image?' value:nearest>, <SpecChoice label:'Interpolate Datacube?' value:linear>, <SpecChoice label:'Render To Use in KML/GeoTIFF' value:NA>, <SpecBool label:'Apply Stretch/Filter?' value:True>, <SpecBool label:'Save Interpolated LCF?' value:False>, <SpecBool label:'Save Mask of Interpolated Pixels?' value:False>]] -> Convert Radiance Cube to Reflectance from Measured Reference Spectrum<SpecGroup>[<SpecSpectrum label:'ROI Spec.' value:spec28>, <SpecBool label:'Measured reflectivity in Percentage?' value:False>, <SpecChoice label:'Scale 100% Reflectivity To:' value:1.0 (floats)>] -> hs_process.crop_single[<SpecPyFloatText label: 'pix_e_ul?' value:342; SpecPyFloatText label: 'pix_n_ul?' value:75; SpecPyFloatText label: 'pix_e_lr?' value:345; SpecPyFloatText label: 'pix_n_lr?' value:78>]
label = Wells_rep2_20180628_16h56m_test_pika_gige_7-Radiance Conversion-Georectify Airborne Datacube-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip
wavelength = {394.6, 396.6528, 398.7056, 400.7584, 402.8112...}
band names = {1, 2, 3, 4, 5...}
```

Notice the "history" tag that describes the history this particular image has gone through to get where it is now. All of the images I work with have this history tag, but if yours do not, then that is probably the cause of the problem.
> Thanks for the information. Can you open up the .hdr file (`E:\hyperspectral_imagery1.hdr`) with a text editor and look at the contents? Here is a sample of the image I used in the example/tutorial:
> 
> ```
> ENVI
> samples = 3
> lines = 3
> bands = 240
> header offset = 0
> file type = ENVI Standard
> data type = 4
> interleave = bip
> byte order = 0
> reflectance scale factor = 1.0
> map info = {UTM, 1.0, 1.0, 441371.00707299996, 4855941.7317699995, 0.04, 0.04, 15, T, WGS-84, units  meters, rotation  0.000}
> history = input cube -> Radiance Conversion<SpecGroup>[<SpecFilename label:'Imager Calibration' value:J:\\\\Shared Work\\\\Data\\\\PikaImagery2_Radiance\\\\Pika II 100106_11 Radiometric Calibration September 2016\\\\8mm\\\\PikaII_100106_11_8mm_RadiometricCalibration_Sept2016.icp>, <SpecCube label:'Dark Noise Cube' value:None>, <SpecBool label:'Auto Remove Dark Noise?' value:True>, <SpecBool label:'Return floating point?' value:False>] -> Georectify Airborne Datacube<SpecGroup>[[<SpecBool label:'Use Flat Earth Altitude?' value:False>, <SpecFloatText label:'Flat Earth Alt (m)' value:340>, <SpecFloatText label:'Field of View (deg)' value:33.0>, <SpecInfo label:'Calculated Resolution (m)' value:Not Available in Batch Mode>, <SpecFloatText label:'Map Resolution (m)' value:0.04>, <SpecBool label:'Select DEM File' value:False>], [<SpecFloatText label:'Sync Offset (s)' value:-0.05>, <SpecFloatText label:'Imager Roll Offset (deg)' value:0>, <SpecFloatText label:'Imager Pitch Offset (deg)' value:2.25>, <SpecFloatText label:'Imager Heading Offset (deg)' value:0>], [<SpecBool label:'Crop Cube to Start/End Lines?' value:False>, <SpecIntText label:'Start Line' value:0>, <SpecIntText label:'End Line' value:1906>], [<SpecBool label:'Select Folder to Save Output' value:False>, <SpecBool label:'Save Products in Source Folder?' value:False>, <SpecBool label:'Generate Full Datacube?' value:True>, <SpecBool label:'Generate KML of Render?' value:True>, <SpecBool label:'Generate GeoTIFF Render?' value:True>, <SpecBool label:'Generate Swath Outline?' value:True>, <SpecChoice label:'Interpolate Image?' value:nearest>, <SpecChoice label:'Interpolate Datacube?' value:linear>, <SpecChoice label:'Render To Use in KML/GeoTIFF' value:NA>, <SpecBool label:'Apply Stretch/Filter?' value:True>, <SpecBool label:'Save Interpolated LCF?' value:False>, <SpecBool label:'Save Mask of Interpolated Pixels?' value:False>]] -> Convert Radiance Cube to Reflectance from Measured Reference Spectrum<SpecGroup>[<SpecSpectrum label:'ROI Spec.' value:spec28>, <SpecBool label:'Measured reflectivity in Percentage?' value:False>, <SpecChoice label:'Scale 100% Reflectivity To:' value:1.0 (floats)>] -> hs_process.crop_single[<SpecPyFloatText label: 'pix_e_ul?' value:342; SpecPyFloatText label: 'pix_n_ul?' value:75; SpecPyFloatText label: 'pix_e_lr?' value:345; SpecPyFloatText label: 'pix_n_lr?' value:78>]
> label = Wells_rep2_20180628_16h56m_test_pika_gige_7-Radiance Conversion-Georectify Airborne Datacube-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip
> wavelength = {394.6, 396.6528, 398.7056, 400.7584, 402.8112...}
> band names = {1, 2, 3, 4, 5...}
> ```
> 
> Notice the "history" tag that describes the history this particular image has gone through to get where it is now. All of the images I work with have this history tag, but if yours do not, then that is probably the cause of the problem.


Yes,  here it is:
```
ENVI
description = {HEADWALL Hyperspec III[HEADWALL Hyperspec III],[HEADWALL Hyperspec III OR]}
samples = 720
lines = 11182
bands = 274
header offset = 0
file type = ENVI Standard
data type = 12
interleave = bsq
sensor type = Unknown
byte order = 0
wavelength units = nm
map info = {Geographic Lat/Lon,1,1,-89.38681637930115,43.29648597737143,3.080735373867722e-7,2.250256219444636e-7,WGS84,units=Degrees}
wavelength = {
399.309
,401.521
,403.732
,405.944
... ...
,998.646
,1000.86
,1003.07
}
;Correct Position = 1
;Firmware = build: Oct 10 2019_09:59:10
;GSDO = 0.025
;HO = 307.0516052246094
;Invert Columns = 1
;Latitude Offset (m) = 0
;Longitude Offset (m) = 0
;Ortho DEM File = C:\Headwall\dem\usgs_ned_13_n44w090_gridfloat.hdr
;Ortho GPS Offset Altitude Offset(m) = 33
;Ortho GPS Offset Front Up Positive = 0
;Ortho GPS Offset Geoid Correction = 1
;Ortho GPS Offset North-East Positive = 1
;Ortho GPS Offset Pitch = -0.4
;Ortho GPS Offset Right Positive = 1
;Ortho GPS Offset Roll = 0
;Ortho GPS Offset Time Offset(ms) = 0
;Ortho GPS Offset Yaw = 3
;Ortho Lems EFL (mm) = 12
;Ortho Lens Pixel Pitch(um) = 7.4
;Ortho Lens Spatial Binning = 1
;Ortho Resolution (m) = 0.025
;Ortho Running Avg. of Altitude = 0
;Ortho Running Avg. of Latitude = 0
;Ortho Running Avg. of Longitude = 0
;Ortho Running Avg. of Pitch = 0
;Ortho Running Avg. of Raw = 0
;Ortho Running Avg. of Roll = 0
;Ortho Running Avg. of X = 0
;Ortho Running Avg. of Y = 0
;Ortho Running Avg. of Yaw = 0
;Ortho Running Avg. of Z = 0
;Ortho Sensor Alpha = 0
;Ortho Sensor Beta = 0
;Ortho Sensor Gamma = 0
;Ortho Sensor Invert Columns = 1
;Row binning = 1
;Serial Number = nHS-213
;averagePixelDisp = 2.21157549723685
;bBinaryMap = 000007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE
;framingCamera = 0
;gain = 3
;hsiii_version = Hyperspec v3.1.0
;invertColumns = 1
;left = 0
;pixel0 = 158.247496218906
;top = 88
```
As I suspected, it looks to be missing the "history" tag. For now until I am able to fix the code, can you manually add ``history = none`` into the .hdr file?

I see all the extra metadata with a leading semicolon - I'm also not sure how ``hs_process`` will behave with these metadata tags. It would be good to get a sample image from you to debug. I've only tested with Resonon imagery.
@Yazhou-Sun, would you mind uploading a single hyperspectral image to [this cloud directory](https://drive.google.com/drive/folders/18AA1HgT1U7CXoSbcAyl7azgoQ3jcryxw?usp=sharing)? I can use it as a test image to make sure the code all works with your images. Thanks!
> @Yazhou-Sun, would you mind uploading a single hyperspectral image to [this cloud directory](https://drive.google.com/drive/folders/18AA1HgT1U7CXoSbcAyl7azgoQ3jcryxw?usp=sharing)? I can use it as a test image to make sure the code all works with your images. Thanks!

Thank you so much for your help! 
I don't have permission to upload files to this drive. would you mind adding me (ysun382@wisc.edu) as an editor or so? 
> > @Yazhou-Sun, would you mind uploading a single hyperspectral image to [this cloud directory](https://drive.google.com/drive/folders/18AA1HgT1U7CXoSbcAyl7azgoQ3jcryxw?usp=sharing)? I can use it as a test image to make sure the code all works with your images. Thanks!
> 
> Thank you so much for your help!
> I don't have permission to upload files to this drive. would you mind adding me ([ysun382@wisc.edu](mailto:ysun382@wisc.edu)) as an editor or so?

Okay, give it a try now - sorry!

Fix patch:
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..68bc525
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,11 @@
+channels:
+    - conda-forge
+dependencies:
+    - python=3.8
+    - geopandas=0.8.*
+    - rasterstats=0.14.*
+    - seaborn=0.11.*
+    - spectral=0.22.*
+    - pip=20.2.*
+    - pip:
+        - -r file:requirements/dev_pip.txt
\ No newline at end of file
diff --git a/hs_process/batch.py b/hs_process/batch.py
index 5bc57bf..75eebda 100644
--- a/hs_process/batch.py
+++ b/hs_process/batch.py
@@ -169,8 +169,8 @@ def _check_processed(self, fname_list, base_dir_out, folder_name,
                 ''.format(len(fname_list_final)))
         if not len(fname_list_final) > 0:
             warnings.warn(msg1, UserWarning, stacklevel=0)
-        else:
-            print(msg2)
+        # else:
+        #     print(msg2)
         time.sleep(0.2)  # when using progress bar, this keeps from splitting lines
         return fname_list_final
 
@@ -877,7 +877,9 @@ def _crop_loop(self, df_plots, gdf, base_dir_out, folder_name,
             fname_hdr = fname + '.hdr'
             self.io.read_cube(fname_hdr, name_long=name_long,
                               name_plot=plot_id_ref, name_short=name_short)
-            self.my_spatial_mod = spatial_mod(self.io.spyfile, gdf)
+            self.my_spatial_mod = spatial_mod(
+                self.io.spyfile, gdf, self.io.base_dir, self.io.name_long,
+                self.io.name_short)
             self.my_spatial_mod.defaults = self.io.defaults
             if base_dir_out is None:
                 dir_out, name_append = self._save_file_setup(
@@ -933,7 +935,9 @@ def _crop_many_read_row(self, row, gdf, method):
         fname_hdr = fname_in + '.hdr'
         self.io.read_cube(fname_hdr, name_long=name_long,
                           name_plot=plot_id_ref, name_short=name_short)
-        self.my_spatial_mod = spatial_mod(self.io.spyfile, gdf)
+        self.my_spatial_mod = spatial_mod(
+            self.io.spyfile, gdf, self.io.base_dir, self.io.name_long,
+            self.io.name_short)
         self.my_spatial_mod.defaults = self.io.defaults
         if method == 'many_gdf':
             df_plots_many = self._many_gdf(cs)
@@ -1035,7 +1039,9 @@ def _execute_crop(self, fname_sheet, fname_list, base_dir_out, folder_name,
                   'pass ``fname_sheet`` with "study" and "date" columns.\n')
             for fname_in in fname_list:
                 self.io.read_cube(fname_in)
-                self.my_spatial_mod = spatial_mod(self.io.spyfile, gdf)
+                self.my_spatial_mod = spatial_mod(
+                    self.io.spyfile, gdf, self.io.base_dir, self.io.name_long,
+                    self.io.name_short))
                 # as long as defaults are set ahead of time, they should carry through
                 # e.g., batch.io.defaults.crop_defaults.n_plots = 40 to limit to 40 plots
                 self.my_spatial_mod.defaults = self.io.defaults
diff --git a/hs_process/spatial_mod.py b/hs_process/spatial_mod.py
index 1c64907..ee74f00 100644
--- a/hs_process/spatial_mod.py
+++ b/hs_process/spatial_mod.py
@@ -4,6 +4,7 @@
 import numpy as np
 import os
 import pandas as pd
+from pyproj import CRS
 from shapely.geometry import Polygon
 import spectral.io.spyfile as SpyFile
 
@@ -16,30 +17,93 @@ class spatial_mod(object):
     Class for manipulating data within the spatial domain
     (e.g., cropping a datacube by a geographical boundary).
     '''
-    def __init__(self, spyfile, gdf=None):
+    def __init__(self, spyfile, gdf=None, **kwargs):
         '''
         spyfile (``SpyFile`` object): The Spectral Python datacube to manipulate.
         gdf (``geopandas.DataFrame``): Polygon data that includes the plot_id and
             its geometry.
+        base_dir (``str``): to be used by the plot gdf attribute data.
+        name_short (``str``): to be used by the plot gdf attribute data.
+        name_long (``str``): to be used by the plot gdf attribute data.
         '''
         self.spyfile = spyfile
         self.gdf = gdf
+        self.base_dir = None
+        self.name_long = None
+        self.name_short = None
+        for k, v in kwargs.items():
+            if k in ['base_dir', 'name_long', 'name_short']:
+                setattr(self, k, v)
 
         self.spy_ps_e = None
         self.spy_ps_n = None
+        self.spy_srid = None
         self.spy_ul_e_srs = None
         self.spy_ul_n_srs = None
         self.tools = None
 
-
         self.defaults = defaults()
-        self.load_spyfile(spyfile)
+        self.load_spyfile(
+            spyfile, base_dir=self.base_dir, name_long=self.name_long,
+            name_short=self.name_short)
+
+    def _get_srid_from_map_info(self, map_info_set):
+        '''
+        Parses the "map info" set of metadata to determine the SRID/EPSG code.
+
+        Images may be in either a geographic or projected coordinate reference
+        system. This function tries to handle any UTM projection, as well as
+        the WGS-84 geographic coordinate system (i.e., latitude/longitude in
+        in units of degrees), however, the projected UTM coordinate reference
+        system is recommended.
 
-    def _create_spyfile_extent_gdf(self, epsg=32615):
+        Note:
+            The "map info" tag of the ENVI header file lists geographic
+            information in the following order:
+                1. Projection name
+                2. Reference (tie point) pixel x location (in file coordinates)
+                3. Reference (tie point) pixel y location (in file coordinates)
+                4. Pixel easting
+                5. Pixel northing
+                6. x pixel size
+                7. y pixel size
+                8. Projection zone (UTM only)
+                9. North or South (UTM only)
+                10. Datum
+                11. Units
+        '''
+        proj_name_list = ['utm', 'geographic', 'wgs']
+        units_list = ['meter', 'feet', 'degree']
+
+        proj_name = [proj_name for proj_name in proj_name_list if proj_name in map_info_set[0].lower()][0]
+        units_set = map_info_set[[i for i, j in enumerate(map_info_set) if 'units' in str(j)][0]]
+        unit = [unit for unit in units_list if unit in units_set.lower()][0]
+
+        if proj_name == 'utm' and unit == 'meter':
+            # Note: EPSG code is retrieved from <map_info_set> index 7 and 8
+            north = True if map_info_set[8] == 'N' or map_info_set[8] == 'T' else False
+            crs_dict = {'proj': proj_name,
+                        'zone': map_info_set[7],
+                        'north': north}
+            crs = CRS.from_dict(crs_dict)
+            srid = int(crs.to_authority()[-1])
+        elif proj_name == 'geographic' and unit == 'degree':
+            srid = 4326
+        else:
+            srid = None
+            print('SRID could not be determined from the "map info" tag. '
+                  'Please make sure the "map info" tag in your .hdr file follows '
+                  'the convention of the ENVI header file format:\n'
+                  'https://www.l3harrisgeospatial.com/docs/enviheaderfiles.html')
+        return srid
+
+    def _create_spyfile_extent_gdf(self, srid=32615):
         '''
+        Creates a geodataframe with a single polygon equal to the extent of the
+        spyfile.
         '''
         metadata = self.spyfile.metadata
-        crs = 'epsg:{0}'.format(epsg)
+        crs = 'epsg:{0}'.format(srid)
         size_x = self.spyfile.shape[1]  # number of pixels
         size_y = self.spyfile.shape[0]
 
@@ -58,7 +122,25 @@ def _create_spyfile_extent_gdf(self, epsg=32615):
         gdf_sp = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])
         return gdf_sp
 
-    def _overlay_gdf(self, gdf, epsg_sp=32615, how='intersection', crop=True):
+    def _check_crs(self, gdf, srid):
+        '''
+        If geopandas CRS is different, reproject.
+
+        NOTE: If CRS is not set but there is also no geometry available, then
+            there won't be changes made.
+        '''
+        if gdf.crs is None and gdf.geometry.isnull().all():
+            return gdf
+        try:
+            if int(gdf.crs.srs.split(':')[-1]) != srid:
+                gdf = gdf.to_crs(epsg=srid)
+        except AttributeError:  # SRS is likely None
+            assert pd.isnull(gdf.crs), 'CRS is not null, yet cannot be determined (?)'
+            gdf.crs = "EPSG:{0}".format(srid)
+            gdf = gdf.to_crs(epsg=srid)
+        return gdf
+
+    def _overlay_gdf(self, gdf, how='intersection', crop=True):
         '''
         Performs a geopandas overlay between the input geodatafram (``gdf``) and
         the extent of ``spyfile``.
@@ -66,7 +148,8 @@ def _overlay_gdf(self, gdf, epsg_sp=32615, how='intersection', crop=True):
         If crop is ``True``, then the polygon bounds will be cropped to the
         extent of the overlay based on the ``how`` parameter.
         '''
-        gdf_sp = self._create_spyfile_extent_gdf(epsg=epsg_sp)
+        gdf_sp = self._create_spyfile_extent_gdf(srid=self.spy_srid)
+        gdf = self._check_crs(gdf, self.spy_srid)
         gdf_filter = gpd.overlay(gdf, gdf_sp, how=how)
         if crop is not True:  # simply return the original plot bounds that intersect
             gdf_filter = gdf[gdf['plot_id'].isin(gdf_filter['plot_id'])].reset_index(drop=True)
@@ -159,9 +242,9 @@ def _find_plots_gdf(self, gdf, plot_id_ref, pix_e_ul, pix_n_ul,
 
             gdf_crop_e_pix = int(abs(plot_srs_e - plot_srs_w) / self.spy_ps_e)
             gdf_crop_n_pix = int(abs(plot_srs_n - plot_srs_s) / self.spy_ps_n)
-            data = [self.tools.base_dir,
-                    self.tools.name_short,
-                    self.tools.name_long,
+            data = [self.base_dir,
+                    self.name_short,
+                    self.name_long,
                     os.path.splitext(self.spyfile.filename)[-1],
                     plot_id, offset_e, offset_n,
                     np.nan, np.nan, np.nan, np.nan,  # buf_X
@@ -352,7 +435,9 @@ def _handle_defaults(self, e_pix, n_pix, e_m, n_m, group='crop',
         if not isinstance(spyfile, SpyFile.SpyFile):
             spyfile = self.spyfile
         else:
-            self.load_spyfile(spyfile)
+            self.load_spyfile(
+                spyfile, base_dir=self.base_dir, name_long=self.name_long,
+                name_short=self.name_short)
 
         if group == 'crop':
             if pd.isnull(e_pix) and pd.isnull(e_m):
@@ -555,7 +640,9 @@ def _crop_many_grid(self, plot_id_ul, pix_e_ul, pix_n_ul,
         if spyfile is None:
             spyfile = self.spyfile
         elif isinstance(spyfile, SpyFile.SpyFile):
-            self.load_spyfile(spyfile)
+            self.load_spyfile(
+                spyfile, base_dir=self.base_dir, name_long=self.name_long,
+                name_short=self.name_short)
 
         msg1 = ('Either crop_size_XX_m or crop_size_XX_pix should be passed. '
                 'Please pass one or the other.')
@@ -804,11 +891,15 @@ def crop_many_gdf(self, spyfile=None, gdf=None, **kwargs):
 
             Read datacube and spatial plot boundaries
 
-            >>> fname_in = r'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip.hdr'
-            >>> fname_gdf = r'F:\\nigo0024\Documents\hs_process_demo\plot_bounds_small\plot_bounds.shp'
+            >>> fname_in = r'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7-Radiance Conversion-Georectify Airborne Datacube-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip.hdr'
+            >>> fname_in = r'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7_nohist-Radiance Conversion-Georectify Airborne Datacube-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip.hdr'
+
+            >>> fname_gdf = r'F:\\nigo0024\Documents\hs_process_demo\plot_bounds.geojson'
             >>> gdf = gpd.read_file(fname_gdf)
             >>> io = hsio(fname_in)
-            >>> my_spatial_mod = spatial_mod(io.spyfile)
+            >>> my_spatial_mod = spatial_mod(
+                    io.spyfile, base_dir=io.base_dir, name_short=io.name_short,
+                    name_long=io.name_long)
             >>> dir_out = os.path.join(io.base_dir, 'spatial_mod', 'crop_many_gdf')
             >>> name_append = '-crop-many-gdf'
 
@@ -888,7 +979,9 @@ def crop_many_gdf(self, spyfile=None, gdf=None, **kwargs):
         kwargs_d = self._check_crop_defaults(**kwargs)
 
         if isinstance(spyfile, SpyFile.SpyFile):
-            self.load_spyfile(spyfile)
+            self.load_spyfile(
+                spyfile, base_dir=self.base_dir, name_long=self.name_long,
+                name_short=self.name_short)
         metadata = self.spyfile.metadata
         if gdf is None:
             gdf = self.gdf
@@ -1029,7 +1122,9 @@ def crop_single(self, pix_e_ul=0, pix_n_ul=0, crop_e_pix=None,
             >>> from hs_process import spatial_mod
             >>> fname_in = r'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip.hdr'
             >>> io = hsio(fname_in)
-            >>> my_spatial_mod = spatial_mod(io.spyfile)
+            >>> my_spatial_mod = spatial_mod(
+                    io.spyfile, base_dir=io.base_dir, name_short=io.name_short,
+                    name_long=io.name_long)
 
             Crop an area with a width (easting) *200 pixels* and a height
             (northing) of *50 pixels*, with a northwest/upper left origin at
@@ -1093,7 +1188,9 @@ def crop_single(self, pix_e_ul=0, pix_n_ul=0, crop_e_pix=None,
             array_crop = spyfile.read_subregion((pix_n_ul, pix_n_lr),
                                                 (pix_e_ul, pix_e_lr))
         elif isinstance(spyfile, SpyFile.SpyFile):
-            self.load_spyfile(spyfile)
+            self.load_spyfile(
+                spyfile, base_dir=self.base_dir, name_long=self.name_long,
+                name_short=self.name_short)
             array_crop = spyfile.read_subregion((pix_n_ul, pix_n_lr),
                                                 (pix_e_ul, pix_e_lr))
         elif isinstance(spyfile, np.ndarray):
@@ -1128,12 +1225,16 @@ def crop_single(self, pix_e_ul=0, pix_n_ul=0, crop_e_pix=None,
         map_info_set = self.tools.modify_meta_set(map_info_set, 4, ul_y_utm)
         metadata['map info'] = map_info_set
 
+        if 'history' not in metadata:  # add history tag to metadata
+            metadata['history'] = '[no prior history]'
+
         hist_str = (" -> hs_process.crop_single[<"
                     "SpecPyFloatText label: 'pix_e_ul?' value:{0}; "
                     "SpecPyFloatText label: 'pix_n_ul?' value:{1}; "
                     "SpecPyFloatText label: 'pix_e_lr?' value:{2}; "
                     "SpecPyFloatText label: 'pix_n_lr?' value:{3}>]"
                     "".format(pix_e_ul, pix_n_ul, pix_e_lr, pix_n_lr))
+
         # If "..crop_single" is already included in the history, remove it
         idx_remove = metadata['history'].find(
                 ' -> hs_process.crop_single[<')
@@ -1155,7 +1256,7 @@ def crop_single(self, pix_e_ul=0, pix_n_ul=0, crop_e_pix=None,
 
         return array_crop, metadata
 
-    def load_spyfile(self, spyfile):
+    def load_spyfile(self, spyfile, **kwargs):
         '''
         Loads a ``SpyFile`` (Spectral Python object) for data access and/or
         manipulation by the ``hstools`` class.
@@ -1163,6 +1264,9 @@ def load_spyfile(self, spyfile):
         Parameters:
             spyfile (``SpyFile`` object): The datacube being accessed and/or
                 manipulated.
+            base_dir (``str``): to be used by the plot gdf attribute data.
+            name_short (``str``): to be used by the plot gdf attribute data.
+            name_long (``str``): to be used by the plot gdf attribute data.
 
         Example:
             Load and initialize the ``hsio`` and ``spatial_mod`` modules
@@ -1171,11 +1275,15 @@ def load_spyfile(self, spyfile):
             >>> from hs_process import spatial_mod
             >>> fname_in = r'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip.hdr'
             >>> io = hsio(fname_in)
-            >>> my_spatial_mod = spatial_mod(io.spyfile)
+            >>> my_spatial_mod = spatial_mod(
+                    io.spyfile, base_dir=io.base_dir, name_short=io.name_short,
+                    name_long=io.name_long)
 
             Load datacube using ``spatial_mod.load_spyfile``
 
-            >>> my_spatial_mod.load_spyfile(io.spyfile)
+            >>> my_spatial_mod.load_spyfile(
+                    io.spyfile, base_dir=io.base_dir, name_short=io.name_short,
+                    name_long=io.name_long)
             >>> my_spatial_mod.spyfile
             Data Source:   'F:\\nigo0024\Documents\hs_process_demo\Wells_rep2_20180628_16h56m_pika_gige_7-Convert Radiance Cube to Reflectance from Measured Reference Spectrum.bip'
         	# Rows:            617
@@ -1185,10 +1293,15 @@ def load_spyfile(self, spyfile):
         	Quantization:  32 bits
         	Data format:   float32
         '''
+        for k, v in kwargs.items():
+            if k in ['base_dir', 'name_long', 'name_short']:
+                setattr(self, k, v)
+
         self.spyfile = spyfile
         self.tools = hstools(spyfile)
         try:
             map_info_set = self.tools.get_meta_set(self.spyfile.metadata['map info'])
+            self.spy_srid = self._get_srid_from_map_info(map_info_set)
             self.spy_ul_e_srs = float(map_info_set[3])
             self.spy_ul_n_srs = float(map_info_set[4])
             self.spy_ps_e = float(map_info_set[5])
@@ -1197,6 +1310,7 @@ def load_spyfile(self, spyfile):
             print('Map information was not able to be loaded from the '
                   '`SpyFile`. Please be sure the metadata contains the "map '
                   'info" tag with accurate geometric information.\n')
+            self.spy_srid = None
             self.spy_ul_e_srs = None
             self.spy_ul_n_srs = None
             self.spy_ps_e = None
diff --git a/hs_process/spec_mod.py b/hs_process/spec_mod.py
index 21c193e..18fe800 100644
--- a/hs_process/spec_mod.py
+++ b/hs_process/spec_mod.py
@@ -58,7 +58,7 @@ def _metadata_derivative(self, metadata, order=1):
         #     metadata_dydx['stdev'] = stdev_dydx
         return metadata_dydx
 
-    def _metadata_mimic(self, sensor, meta_bands_sensor):
+    def _metadata_mimic(self, sensor, meta_bands_sensor, fwhm_dict):
         '''Modifies metadata for spectral_mimic() function.'''
         metadata = deepcopy(self.tools.spyfile.metadata.copy())
         band_names = list(meta_bands_sensor.keys())
@@ -71,8 +71,10 @@ def _metadata_mimic(self, sensor, meta_bands_sensor):
         wavelength = list(meta_bands_sensor.values())
         band_str = '{' + ', '.join(str(b) for b in band_names) + '}'
         wavelength_str = '{' + ', '.join(str(wl) for wl in wavelength) + '}'
+        fwhm_str = '{' + ', '.join(str(fwhm_dict[b]) for b in band_names) + '}'
         metadata['band names'] = band_str
         metadata['wavelength'] = wavelength_str
+        metadata['fwhm'] = fwhm_str
         # self.tools.spyfile.metadata = metadata
         return metadata
 
@@ -146,6 +148,22 @@ def _mimic_get_response(self, sensor, df_band_response, col_wl):
         band_names.remove('wl_nm')
         return df, band_names
 
+    def _mimic_get_fwhm(self, df):
+        '''
+        Gets the FWHM off the band responses in ``df``.
+
+        Parameters:
+            df: The band response dataframe (first column should be "wl_nm"
+                indicating the wavelength).
+        '''
+        fwhm_dict = {}
+        for c in df.columns:
+            if c == 'wl_nm':
+                continue
+            df_half_max = df[df[c] >= 0.5]
+            fwhm_dict[c] = df_half_max['wl_nm'].max() - df_half_max['wl_nm'].min()
+        return fwhm_dict
+
     def _resample_get_bin_indices(self, bandwidth, bins_n):
         '''
         Gets an array of indices that indicate which bands should be averaged/
@@ -524,6 +542,7 @@ def spectral_mimic(self, sensor='sentinel-2a', df_band_response=None,
             entire image).
 
             >>> import seaborn as sns
+            >>> from ast import literal_eval
             >>> spy_hs = my_spec_mod.spyfile.open_memmap()  # datacube before smoothing
             >>> meta_bands = list(io.tools.meta_bands.values())
             >>> meta_bands_s2a = sorted([float(i) for i in literal_eval(metadata_s2a['wavelength'])])
@@ -612,7 +631,8 @@ def spectral_mimic(self, sensor='sentinel-2a', df_band_response=None,
                 bands_remove[idx] = band_name
         for idx, band_name in bands_remove.items():
             array_multi = np.delete(array_multi, idx, axis=2)
-        metadata = self._metadata_mimic(sensor, meta_bands_sensor)
+        fwhm_dict = self._mimic_get_fwhm(df)
+        metadata = self._metadata_mimic(sensor, meta_bands_sensor, fwhm_dict)
         return array_multi, metadata
 
     def spectral_resample(self, bandwidth=None, bins_n=None, spyfile=None):
diff --git a/hs_process/utilities.py b/hs_process/utilities.py
index cf9a05d..b5b92a2 100644
--- a/hs_process/utilities.py
+++ b/hs_process/utilities.py
@@ -551,7 +551,7 @@ def _parse_fname(self, fname_hdr=None, str_plot='plot_', overwrite=True,
 
         if name_plot is not None:
             self.name_plot = name_plot
-        elif overwrite is False and self.name_plot is None:
+        elif overwrite is False and self.name_plot is None:  # don't have name_plot yet, so see if we can find it
             self.name_plot = self._parse_fname_plot(str_plot)
         else:
             pass
@@ -1544,18 +1544,20 @@ class hstools(object):
         spyfile (``SpyFile`` object): The datacube being accessed and/or
             manipulated.
     '''
+    # def __init__(self, spyfile):
     def __init__(self, spyfile):
         msg = ('Pleae load a SpyFile (Spectral Python object)')
         assert spyfile is not None, msg
 
         self.spyfile = spyfile
 
-        self.fname_in = None
-        self.fname_hdr = None
-        self.base_name = None
-        self.name_short = None
-        self.name_long = None
-        self.name_plot = None
+        # self.fname_in = None
+        # self.fname_hdr = None
+        # self.base_dir = None
+        # self.base_name = None
+        # self.name_short = None
+        # self.name_long = None
+        # self.name_plot = None
         self.meta_bands = None
 
         self.load_spyfile(spyfile)
@@ -1598,107 +1600,107 @@ def _get_meta_bands(self, spyfile=None, metadata=None):
                     meta_bands[band_names[idx]] = float(wl_names[idx])
         self.meta_bands = meta_bands
 
-    def _parse_fname_plot(self, str_plot):
-        '''
-        Code for parsing ``name_plot`` (numeric text following ``str_plot``).
-        '''
-        s = self.name_short
-        if str_plot in s and '_pika' in s:
-            name_plot = s[s.find(str_plot) + len(str_plot):s.find('_pika')]
-        elif str_plot in s and '_pika' not in s:
-            name_plot = s[s.find(str_plot) + len(str_plot):s.find('-')]
-        else:
-            name_plot = self.name_short.rsplit('_', 1)[1]
-
-        if len(name_plot) > 12:  # then it must have gone wrong
-            name_plot = self.name_short.rsplit('_', 1)[1]
-        if len(name_plot) == 0:  # then '_pika' must not
-            name_plot = self.name_short.rsplit('_', 1)[1]
-
-        try:
-            int(name_plot)
-        except ValueError:  # give up..
-            msg = ('Cannot determine the plot name from the image filename. '
-                   'Setting `hsio.name_plot` to `None`. If this image is for '
-                   'a particular plot, please set `hsio.name_plot; otherwise, '
-                   'ignore this warning.\n')
-            warnings.warn(msg, UserWarning)
-
-            name_plot = None
-        return name_plot
-
-    def _parse_fname(self, fname_hdr=None, str_plot='plot_', overwrite=True,
-                     name_long=None, name_short=None, name_plot=None):
-        '''
-        Parses the filename for ``name_long`` (text after the first dash,
-        inclusive), ``name_short`` (text before the first dash), and
-        ``name_plot`` (numeric text following ``str_plot``).
-
-        Parameters:
-            fname_hdr (``str``): input filename.
-            str_plot (``str``): text to search for that precedes the numeric
-                text that describes the plot number.
-            overwrite (``bool``): whether the class instances of ``name_long``,
-                ``name_short``, and ``name_plot`` should be overwritten based
-                on ``fname_in`` (default: ``True``).
-        '''
-        if fname_hdr is None:
-            fname_hdr = self.spyfile.filename + '.hdr'
-        if os.path.splitext(fname_hdr)[1] == '.hdr':  # modify self.fname_in based on new file
-            fname_in = os.path.splitext(fname_hdr)[0]
-        else:
-            fname_hdr = fname_hdr + '.hdr'
-            fname_in = os.path.splitext(fname_hdr)[0]
-        self.fname_in = fname_in
-        self.fname_hdr = fname_hdr
-
-        self.base_dir = os.path.dirname(fname_in)
-#        base_name = os.path.basename(fname_in)
-        base_name = os.path.basename(os.path.splitext(fname_in)[0])
-        self.base_name = base_name
-
-        if overwrite is True:
-            if '-' in base_name:
-                self.name_long = base_name[base_name.find(
-                        '-', base_name.rfind('_')):]
-                self.name_short = base_name[:base_name.find(
-                        '-', base_name.rfind('_'))]
-            else:
-                # if name_long does not have ext, it can be just blank
-                self.name_long = ''
-                # and name_short can be base_name
-                self.name_short = base_name
-            self.name_plot = self._parse_fname_plot(str_plot)
-
-        if name_long is not None:
-            self.name_long = name_long
-        elif overwrite is False and self.name_long is None:
-            if '-' in base_name:
-                self.name_long = base_name[base_name.find(
-                        '-', base_name.rfind('_')):]
-            else:
-                # if name_long does not have ext, it can be just blank
-                self.name_long = ''
-        else:
-            pass
-
-        if name_short is not None:
-            self.name_short = name_short
-        elif overwrite is False and self.name_short is None:
-            if '-' in base_name:
-                self.name_short = base_name[:base_name.find(
-                        '-', base_name.rfind('_'))]
-            else:
-                self.name_short = base_name
-        else:
-            pass
-
-        if name_plot is not None:
-            self.name_plot = name_plot
-        elif overwrite is False and self.name_plot is None:
-            self.name_plot = self._parse_fname_plot(str_plot)
-        else:
-            pass
+#     def _parse_fname_plot(self, str_plot):
+#         '''
+#         Code for parsing ``name_plot`` (numeric text following ``str_plot``).
+#         '''
+#         s = self.name_short
+#         if str_plot in s and '_pika' in s:
+#             name_plot = s[s.find(str_plot) + len(str_plot):s.find('_pika')]
+#         elif str_plot in s and '_pika' not in s:
+#             name_plot = s[s.find(str_plot) + len(str_plot):s.find('-')]
+#         else:
+#             name_plot = self.name_short.rsplit('_', 1)[1]
+
+#         if len(name_plot) > 12:  # then it must have gone wrong
+#             name_plot = self.name_short.rsplit('_', 1)[1]
+#         if len(name_plot) == 0:  # then '_pika' must not
+#             name_plot = self.name_short.rsplit('_', 1)[1]
+
+#         try:
+#             int(name_plot)
+#         except ValueError:  # give up..
+#             msg = ('Cannot determine the plot name from the image filename. '
+#                    'Setting `hsio.name_plot` to `None`. If this image is for '
+#                    'a particular plot, please set `hsio.name_plot; otherwise, '
+#                    'ignore this warning.\n')
+#             warnings.warn(msg, UserWarning)
+
+#             name_plot = None
+#         return name_plot
+
+#     def _parse_fname(self, fname_hdr=None, str_plot='plot_', overwrite=True,
+#                      name_long=None, name_short=None, name_plot=None):
+#         '''
+#         Parses the filename for ``name_long`` (text after the first dash,
+#         inclusive), ``name_short`` (text before the first dash), and
+#         ``name_plot`` (numeric text following ``str_plot``).
+
+#         Parameters:
+#             fname_hdr (``str``): input filename.
+#             str_plot (``str``): text to search for that precedes the numeric
+#                 text that describes the plot number.
+#             overwrite (``bool``): whether the class instances of ``name_long``,
+#                 ``name_short``, and ``name_plot`` should be overwritten based
+#                 on ``fname_in`` (default: ``True``).
+#         '''
+#         if fname_hdr is None:
+#             fname_hdr = self.spyfile.filename + '.hdr'
+#         if os.path.splitext(fname_hdr)[1] == '.hdr':  # modify self.fname_in based on new file
+#             fname_in = os.path.splitext(fname_hdr)[0]
+#         else:
+#             fname_hdr = fname_hdr + '.hdr'
+#             fname_in = os.path.splitext(fname_hdr)[0]
+#         self.fname_in = fname_in
+#         self.fname_hdr = fname_hdr
+
+#         self.base_dir = os.path.dirname(fname_in)
+# #        base_name = os.path.basename(fname_in)
+#         base_name = os.path.basename(os.path.splitext(fname_in)[0])
+#         self.base_name = base_name
+
+#         if overwrite is True:
+#             if '-' in base_name:
+#                 self.name_long = base_name[base_name.find(
+#                         '-', base_name.rfind('_')):]
+#                 self.name_short = base_name[:base_name.find(
+#                         '-', base_name.rfind('_'))]
+#             else:
+#                 # if name_long does not have ext, it can be just blank
+#                 self.name_long = ''
+#                 # and name_short can be base_name
+#                 self.name_short = base_name
+#             self.name_plot = self._parse_fname_plot(str_plot)
+
+#         if name_long is not None:
+#             self.name_long = name_long
+#         elif overwrite is False and self.name_long is None:
+#             if '-' in base_name:
+#                 self.name_long = base_name[base_name.find(
+#                         '-', base_name.rfind('_')):]
+#             else:
+#                 # if name_long does not have ext, it can be just blank
+#                 self.name_long = ''
+#         else:
+#             pass
+
+#         if name_short is not None:
+#             self.name_short = name_short
+#         elif overwrite is False and self.name_short is None:
+#             if '-' in base_name:
+#                 self.name_short = base_name[:base_name.find(
+#                         '-', base_name.rfind('_'))]
+#             else:
+#                 self.name_short = base_name
+#         else:
+#             pass
+
+#         if name_plot is not None:
+#             self.name_plot = name_plot
+#         elif overwrite is False and self.name_plot is None:  # don't have name_plot yet, so see if we can find it
+#             self.name_plot = self._parse_fname_plot(str_plot)
+#         else:
+#             pass
 
     def clean_md_sets(self, metadata=None):
         '''
@@ -2392,8 +2394,9 @@ def load_spyfile(self, spyfile):
         '''
         self.spyfile = spyfile
         self._get_meta_bands(spyfile)
-        self._parse_fname(fname_hdr=None, str_plot='plot_', overwrite=True,
-                          name_long=None, name_short=None, name_plot=None)
+        # self._parse_fname(fname_hdr=None, str_plot='plot_', overwrite=True,
+        #                   name_long=None, name_short=None, name_plot=None)
+        self.base_dir = os.path.dirname(self.spyfile.filename)
 
     def mask_array(self, array, metadata, thresh=None, percentile=None,
                    side='lower'):
diff --git a/requirements/dev_pip.txt b/requirements/dev_pip.txt
new file mode 100644
index 0000000..b2b248a
--- /dev/null
+++ b/requirements/dev_pip.txt
@@ -0,0 +1,1 @@
+git+https://github.com/maroba/findiff.git
\ No newline at end of file


Test patch:
diff --git a/hs_process/tests/test_hsio.py b/hs_process/tests/test_hsio.py
index 0127220..79fadc0 100644
--- a/hs_process/tests/test_hsio.py
+++ b/hs_process/tests/test_hsio.py
@@ -232,8 +232,10 @@ def setUp(self):
         self.fname_hdr = FILENAME_HDR
         self.test_dir = tempfile.mkdtemp()
         self.io = hsio(FILENAME_HDR)
-        self.my_spatial_mod = spatial_mod(self.io.spyfile)
-
+        self.my_spatial_mod = spatial_mod(
+            self.io.spyfile, base_dir=self.io.base_dir,
+            name_short=self.io.name_short,
+            name_long=self.io.name_long)
     def tearDown(self):
         '''
         This tearDown function will be called after each test method is run
@@ -292,7 +294,10 @@ def setUp(self):
         self.test_dir = tempfile.mkdtemp()
         self.fname_hdr = FILENAME_HDR
         self.io = hsio(FILENAME_HDR)
-        self.my_spatial_mod = spatial_mod(self.io.spyfile)
+        self.my_spatial_mod = spatial_mod(
+            self.io.spyfile, base_dir=self.io.base_dir,
+            name_short=self.io.name_short,
+            name_long=self.io.name_long)
 
     def tearDown(self):
         '''