diff --git a/other/datasets/placenames/oqaasileriffik_placenames.cpg b/other/datasets/placenames/oqaasileriffik_placenames.cpg
deleted file mode 100644
index 3ad133c..0000000
--- a/other/datasets/placenames/oqaasileriffik_placenames.cpg
+++ /dev/null
@@ -1 +0,0 @@
-UTF-8
\ No newline at end of file
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.dbf b/other/datasets/placenames/oqaasileriffik_placenames.dbf
deleted file mode 100644
index f36df80..0000000
Binary files a/other/datasets/placenames/oqaasileriffik_placenames.dbf and /dev/null differ
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.prj b/other/datasets/placenames/oqaasileriffik_placenames.prj
deleted file mode 100644
index 212c2e0..0000000
--- a/other/datasets/placenames/oqaasileriffik_placenames.prj
+++ /dev/null
@@ -1 +0,0 @@
-PROJCS["WGS_1984_UTM_Zone_24N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-39.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.sbn b/other/datasets/placenames/oqaasileriffik_placenames.sbn
deleted file mode 100644
index aeb00ec..0000000
Binary files a/other/datasets/placenames/oqaasileriffik_placenames.sbn and /dev/null differ
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.sbx b/other/datasets/placenames/oqaasileriffik_placenames.sbx
deleted file mode 100644
index ec6d89c..0000000
Binary files a/other/datasets/placenames/oqaasileriffik_placenames.sbx and /dev/null differ
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.shp b/other/datasets/placenames/oqaasileriffik_placenames.shp
deleted file mode 100644
index afc7374..0000000
Binary files a/other/datasets/placenames/oqaasileriffik_placenames.shp and /dev/null differ
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.shp.xml b/other/datasets/placenames/oqaasileriffik_placenames.shp.xml
deleted file mode 100644
index 288bb51..0000000
--- a/other/datasets/placenames/oqaasileriffik_placenames.shp.xml
+++ /dev/null
@@ -1,2 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<metadata xml:lang="en"><Esri><CreaDate>20191104</CreaDate><CreaTime>12103400</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>FALSE</SyncOnce><DataProperties><itemProps><itemName Sync="TRUE">lakenames_pts</itemName><imsContentType Sync="TRUE">002</imsContentType><itemSize Sync="TRUE">0.000</itemSize><itemLocation><linkage Sync="TRUE">file://\\ASIAQWS012\D$\python_workspace\ice_marginal_lakes\lakenames\lakenames_pts.shp</linkage><protocol Sync="TRUE">Local Area Network</protocol></itemLocation></itemProps><coordRef><type Sync="TRUE">Projected</type><geogcsn Sync="TRUE">GCS_WGS_1984</geogcsn><csUnits Sync="TRUE">Linear Unit: Meter (1.000000)</csUnits><projcsn Sync="TRUE">WGS_1984_UTM_Zone_24N</projcsn><peXml Sync="TRUE">&lt;ProjectedCoordinateSystem xsi:type='typens:ProjectedCoordinateSystem' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/2.1.0'&gt;&lt;WKT&gt;PROJCS[&amp;quot;WGS_1984_UTM_Zone_24N&amp;quot;,GEOGCS[&amp;quot;GCS_WGS_1984&amp;quot;,DATUM[&amp;quot;D_WGS_1984&amp;quot;,SPHEROID[&amp;quot;WGS_1984&amp;quot;,6378137.0,298.257223563]],PRIMEM[&amp;quot;Greenwich&amp;quot;,0.0],UNIT[&amp;quot;Degree&amp;quot;,0.0174532925199433]],PROJECTION[&amp;quot;Transverse_Mercator&amp;quot;],PARAMETER[&amp;quot;False_Easting&amp;quot;,500000.0],PARAMETER[&amp;quot;False_Northing&amp;quot;,0.0],PARAMETER[&amp;quot;Central_Meridian&amp;quot;,-39.0],PARAMETER[&amp;quot;Scale_Factor&amp;quot;,0.9996],PARAMETER[&amp;quot;Latitude_Of_Origin&amp;quot;,0.0],UNIT[&amp;quot;Meter&amp;quot;,1.0],AUTHORITY[&amp;quot;EPSG&amp;quot;,32624]]&lt;/WKT&gt;&lt;XOrigin&gt;-5120900&lt;/XOrigin&gt;&lt;YOrigin&gt;-9998100&lt;/YOrigin&gt;&lt;XYScale&gt;450445547.3910538&lt;/XYScale&gt;&lt;ZOrigin&gt;-100000&lt;/ZOrigin&gt;&lt;ZScale&gt;10000&lt;/ZScale&gt;&lt;MOrigin&gt;-100000&lt;/MOrigin&gt;&lt;MScale&gt;10000&lt;/MScale&gt;&lt;XYTolerance&gt;0.001&lt;/XYTolerance&gt;&lt;ZTolerance&gt;0.001&lt;/ZTolerance&gt;&lt;MTolerance&gt;0.001&lt;/MTolerance&gt;&lt;HighPrecision&gt;true&lt;/HighPrecision&gt;&lt;WKID&gt;32624&lt;/WKID&gt;&lt;LatestWKID&gt;32624&lt;/LatestWKID&gt;&lt;/ProjectedCoordinateSystem&gt;</peXml></coordRef><lineage><Process ToolSource="c:\program files\arcgis\pro\Resources\ArcToolbox\toolboxes\Data Management Tools.tbx\CopyFeatures" Date="20191104" Time="121049">CopyFeatures 20190218_Stednavne_point D:\python_workspace\ice_marginal_lakes\lakenames\lakenames_pts.shp # # # #</Process></lineage></DataProperties><SyncDate>20191104</SyncDate><SyncTime>12103400</SyncTime><ModDate>20191104</ModDate><ModTime>12103400</ModTime></Esri><dataIdInfo><envirDesc Sync="TRUE">Microsoft Windows 10 Version 10.0 (Build 18362) ; Esri ArcGIS 12.1.1.10257</envirDesc><dataLang><languageCode value="eng" Sync="TRUE"></languageCode><countryCode value="GBR" Sync="TRUE"></countryCode></dataLang><idCitation><resTitle Sync="TRUE">lakenames_pts</resTitle><presForm><PresFormCd value="005" Sync="TRUE"></PresFormCd></presForm></idCitation><spatRpType><SpatRepTypCd value="001" Sync="TRUE"></SpatRepTypCd></spatRpType></dataIdInfo><mdLang><languageCode value="eng" Sync="TRUE"></languageCode><countryCode value="GBR" Sync="TRUE"></countryCode></mdLang><mdChar><CharSetCd value="004" Sync="TRUE"></CharSetCd></mdChar><distInfo><distFormat><formatName Sync="TRUE">Shapefile</formatName></distFormat><distTranOps><transSize Sync="TRUE">0.000</transSize></distTranOps></distInfo><mdHrLv><ScopeCd value="005" Sync="TRUE"></ScopeCd></mdHrLv><mdHrLvName Sync="TRUE">dataset</mdHrLvName><refSysInfo><RefSystem><refSysID><identCode code="32624" Sync="TRUE"></identCode><idCodeSpace Sync="TRUE">EPSG</idCodeSpace><idVersion Sync="TRUE">2.1(3.0.1)</idVersion></refSysID></RefSystem></refSysInfo><spatRepInfo><VectSpatRep><geometObjs Name="lakenames_pts"><geoObjTyp><GeoObjTypCd value="004" Sync="TRUE"></GeoObjTypCd></geoObjTyp><geoObjCnt Sync="TRUE">0</geoObjCnt></geometObjs><topLvl><TopoLevCd value="001" Sync="TRUE"></TopoLevCd></topLvl></VectSpatRep></spatRepInfo><spdoinfo><ptvctinf><esriterm Name="lakenames_pts"><efeatyp Sync="TRUE">Simple</efeatyp><efeageom code="1" Sync="TRUE"></efeageom><esritopo Sync="TRUE">FALSE</esritopo><efeacnt Sync="TRUE">0</efeacnt><spindex Sync="TRUE">FALSE</spindex><linrefer Sync="TRUE">FALSE</linrefer></esriterm></ptvctinf></spdoinfo><eainfo><detailed Name="lakenames_pts"><enttyp><enttypl Sync="TRUE">lakenames_pts</enttypl><enttypt Sync="TRUE">Feature Class</enttypt><enttypc Sync="TRUE">0</enttypc></enttyp><attr><attrlabl Sync="TRUE">FID</attrlabl><attalias Sync="TRUE">FID</attalias><attrtype Sync="TRUE">OID</attrtype><attwidth Sync="TRUE">4</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale><attrdef Sync="TRUE">Internal feature number.</attrdef><attrdefs Sync="TRUE">Esri</attrdefs><attrdomv><udom Sync="TRUE">Sequential unique whole numbers that are automatically generated.</udom></attrdomv></attr><attr><attrlabl Sync="TRUE">Shape</attrlabl><attalias Sync="TRUE">Shape</attalias><attrtype Sync="TRUE">Geometry</attrtype><attwidth Sync="TRUE">0</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale><attrdef Sync="TRUE">Feature geometry.</attrdef><attrdefs Sync="TRUE">Esri</attrdefs><attrdomv><udom Sync="TRUE">Coordinates defining the features.</udom></attrdomv></attr><attr><attrlabl Sync="TRUE">ID</attrlabl><attalias Sync="TRUE">ID</attalias><attrtype Sync="TRUE">Integer</attrtype><attwidth Sync="TRUE">10</attwidth><atprecis Sync="TRUE">10</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Ny_grønla</attrlabl><attalias Sync="TRUE">Ny_grønla</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">150</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Gammel_gr</attrlabl><attalias Sync="TRUE">Gammel_gr</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">150</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Dansk</attrlabl><attalias Sync="TRUE">Dansk</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">150</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Alternativ</attrlabl><attalias Sync="TRUE">Alternativ</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">150</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Genstand</attrlabl><attalias Sync="TRUE">Genstand</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">30</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Kilde</attrlabl><attalias Sync="TRUE">Kilde</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Kommune</attrlabl><attalias Sync="TRUE">Kommune</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">20</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Noter</attrlabl><attalias Sync="TRUE">Noter</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Oprettelse</attrlabl><attalias Sync="TRUE">Oprettelse</attalias><attrtype Sync="TRUE">Date</attrtype><attwidth Sync="TRUE">8</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Senest_ret</attrlabl><attalias Sync="TRUE">Senest_ret</attalias><attrtype Sync="TRUE">Date</attrtype><attwidth Sync="TRUE">8</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Autorisere</attrlabl><attalias Sync="TRUE">Autorisere</attalias><attrtype Sync="TRUE">Integer</attrtype><attwidth Sync="TRUE">5</attwidth><atprecis Sync="TRUE">5</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Autorise00</attrlabl><attalias Sync="TRUE">Autorise00</attalias><attrtype Sync="TRUE">Date</attrtype><attwidth Sync="TRUE">8</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Projektsty</attrlabl><attalias Sync="TRUE">Projektsty</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">100</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">KMS_id</attrlabl><attalias Sync="TRUE">KMS_id</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">5</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Bredde</attrlabl><attalias Sync="TRUE">Bredde</attalias><attrtype Sync="TRUE">Double</attrtype><attwidth Sync="TRUE">19</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Længde</attrlabl><attalias Sync="TRUE">Længde</attalias><attrtype Sync="TRUE">Double</attrtype><attwidth Sync="TRUE">19</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Blad</attrlabl><attalias Sync="TRUE">Blad</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">5</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">DMS_længd</attrlabl><attalias Sync="TRUE">DMS_længd</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">10</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">DMS_bredde</attrlabl><attalias Sync="TRUE">DMS_bredde</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">10</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">KMS_kartog</attrlabl><attalias Sync="TRUE">KMS_kartog</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">15</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Navngiver</attrlabl><attalias Sync="TRUE">Navngiver</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Navngivet</attrlabl><attalias Sync="TRUE">Navngivet</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Betydning</attrlabl><attalias Sync="TRUE">Betydning</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr><attr><attrlabl Sync="TRUE">Historie</attrlabl><attalias Sync="TRUE">Historie</attalias><attrtype Sync="TRUE">String</attrtype><attwidth Sync="TRUE">254</attwidth><atprecis Sync="TRUE">0</atprecis><attscale Sync="TRUE">0</attscale></attr></detailed></eainfo><mdDateSt Sync="TRUE">20191104</mdDateSt></metadata>
diff --git a/other/datasets/placenames/oqaasileriffik_placenames.shx b/other/datasets/placenames/oqaasileriffik_placenames.shx
deleted file mode 100644
index 8351348..0000000
Binary files a/other/datasets/placenames/oqaasileriffik_placenames.shx and /dev/null differ
diff --git a/src/griml/stats/general_stats.py b/src/griml/stats/general_stats.py
index dc457be..69896b3 100644
--- a/src/griml/stats/general_stats.py
+++ b/src/griml/stats/general_stats.py
@@ -2,9 +2,11 @@
 
 import numpy as np
 import geopandas as gp
+import glob
+from pathlib import Path
 
 
-def general_stats(geofile, outfile): 
+def general_stats(file1, outfile=None): 
     '''Calculate general statistics on a lake inventory file
     
     Parameters
@@ -14,81 +16,98 @@ def general_stats(geofile, outfile):
     outfile : str
         Outputted file name for general statistics
     '''
-    
-    #Read file as geopandas index
+    print('Retrieving lake data...') 
     geofile = gp.read_file(file1)
-    agg_geofile = geofile.dissolve(by='unique_id')
+    geofile = geofile.sort_values('region')    
+    uniquebasin = list(geofile['region'].unique())
+    uniquesource = list(geofile['source'].unique())
+    
+    print('Aggregating lake data...')
+    agg_geofile = geofile.dissolve(by='lake_id')
+    agg_geofile['area_sqkm']=[g.area/10**6 for g in list(agg_geofile['geometry'])]
+    agg_geofile['length_km']=[g.length/1000 for g in list(agg_geofile['geometry'])]
+    agg_geofile.sort_values('region')  
 
-    print('\nRetrieving lake data...')
-    
-    #Get data from columns
-    geofile = geofile.sort_values('drainageba')
-    geofile_id = geofile['unique_id'].tolist()              #Get lake IDs
-    geofile_uncertain = geofile['certainty'].tolist()       #Get lake uncertainty
-    geofile_source = geofile['source'].tolist()             #Get lake source
-    geofile_basin = geofile['drainageba'].tolist()          #Get lake location
-    geofile['area'] = geofile['geometry'].area/10**6
-    
-    #Get unique values for basin area and satellite source 
-    uniquebasin = list(set(geofile_basin))
-    uniquesource = list(set(geofile_source))
-    
-    print('\nAggregating lake data...')
-    
-    #Get data from columns
-    agg_geofile['area'] = geofile['geometry'].area
-    agg_geofile['length'] = geofile['geometry'].length
-    agg_geofile.sort_values('drainageba')  
-    aggfile_basin = agg_geofile['drainageba'].tolist()    #Get lake location
-    aggfile_area = agg_geofile['area'].tolist()           #Get lake area
-    aggfile_areakm = []
-    for i in aggfile_area:
-        aggfile_areakm.append(i/10**6)
-    
-    label=['CW', 'CE', 'IC', 'NE', 'NO', 'NW', 'SE', 'SW']
-               
-    print('\nWriting general stats file...')
-            
-    #Open stats file    
-    f = open(outtxt1, 'w+')
-    
-    #Total number of lakes in dataset and total number of unique lakes
-    f.write('Total number of detected lakes , ' + str(len(geofile_id)) + '\n')
-    f.write('Total number of unique lakes , ' + str(max(geofile_id)) + '\n\n')
-    
-    #Total lake count for each sector
-    for i in range(len(label)):
-        f.write(label[i] + ' total lake count , ' + str(geofile_basin.count(label[i])) + '\n')
-    f.write('\n')
-    
-    #Unique lake count for each sector
-    for i in range(len(label)):
-        f.write(label[i] + ' unique lake count , ' + str(aggfile_basin.count(label[i])) + '\n')
-    f.write('\n')
-    
-    #Source count
+    print('\nStarting general statistics...')    
+    readout=''
+    readout=readout+'Total number of detected lakes: ' + str(len(list(geofile['lake_id']))) + '\n'
+    readout=readout+'Total number of unique lakes: ' + str(geofile['lake_id'].max()) + '\n\n'
+
+    # Total lake count for each ice sheet sector
+    geofile_icesheet = geofile[geofile['margin'] == 'ICE_SHEET']
+    agggeofile_icesheet = agg_geofile[agg_geofile['margin'] == 'ICE_SHEET']
+    readout=readout+'Ice Sheet total lake count\n' 
+    for i in uniquebasin:
+        readout=readout+i + ' total lake count: ' + str(geofile_icesheet['region'].value_counts()[i]) + '\n'
+    readout=readout+'Ice sheet total lake count: ' + str(len(geofile_icesheet))
+    readout=readout+'\n\n'
+    
+    # Unique lake count for each ice sheet sector
+    readout=readout+'Ice Sheet unique lake count\n' 
+    for i in uniquebasin:
+        readout=readout+i + ' unique lake count: ' + str(agggeofile_icesheet['region'].value_counts()[i]) + '\n'
+    readout=readout+'Ice sheet total unique lake count: ' + str(len(agggeofile_icesheet ))
+    readout=readout+'\n\n'
+
+    # Total lake count for each ice cap sector
+    geofile_icecap = geofile[geofile['margin'] == 'ICE_CAP']
+    agggeofile_icecap = agg_geofile[agg_geofile['margin'] == 'ICE_CAP']
+    readout=readout+'Ice Cap total lake count\n'
+    for i in uniquebasin:
+        readout=readout+i + ' total lake count: ' + str(geofile_icecap['region'].value_counts()[i]) + '\n'
+    readout=readout+'Ice cap total lake count: ' + str(len(geofile_icecap))
+    readout=readout+'\n\n'
+    
+    # Unique lake count for each ice cap sector
+    readout=readout+'Ice Cap unique lake count\n' 
+    for i in uniquebasin:
+        readout=readout+i+ ' unique lake count: ' + str(agggeofile_icecap['region'].value_counts()[i]) + '\n'
+    readout=readout+'Ice cap total unique lake count: ' + str(len(agggeofile_icecap))
+    readout=readout+'\n\n'
+
+    # Source count
     for i in uniquesource:
-        f.write('Lakes detected using ' + i + ',' + str(geofile_source.count(i)) + '\n')
-    f.write('\n')
-    
-    #Min, max and average lake area
-    f.write('Min. lake area (km), ' + str(min(aggfile_areakm)) + '\n')
-    f.write('Max. lake area (km), ' + str(max(aggfile_areakm)) + '\n')
-    f.write('Average lake area (km), ' + str(np.average(aggfile_areakm)) + '\n')
-    f.write('\n')
-    
-    #Min, max and average uncertainty
-    f.write('Min. uncertainty , ' + str(min(geofile_uncertain)) + '\n')
-    f.write('Max. uncertainty , ' + str(max(geofile_uncertain)) + '\n')
-    f.write('Average uncertainty , ' + str(np.average(geofile_uncertain)) + '\n\n')
-    
-    f.close()
+        readout=readout+'Lakes detected using ' + i + ': ' + str(geofile['source'].value_counts()[i]) + '\n'
+    readout=readout+'\n'
+    
+    # Min, max and average lake area
+    readout=readout+'Min. lake area (km): ' + str(agg_geofile['area_sqkm'].min()) + '\n'
+    readout=readout+'Max. lake area (km): ' + str(agg_geofile['area_sqkm'].max()) + '\n'
+    readout=readout+'Average lake area (km): ' + str(agg_geofile['area_sqkm'].max()) + '\n\n'
+ 
+    # Min, max and average uncertainty
+    readout=readout+'Min. uncertainty: ' + str(agg_geofile['certainty'].min()) + '\n'
+    readout=readout+'Max. uncertainty: ' + str(agg_geofile['certainty'].max()) + '\n'
+    readout=readout+'Average uncertainty: ' + str(agg_geofile['certainty'].mean()) + '\n\n'
+
+    print(readout)
+                  
+    # Open stats file   
+    if outfile is not None:
+        print('Writing general stats file...')
+        f = open(outfile, 'w+')
+        f.write(readout)    
+        f.close()
+    
+    return readout
 
 if __name__ == "__main__": 
     
-    #Define input file and location 
-    workspace1 = '/home/pho/python_workspace/GrIML/other/'
-    file1 = workspace1 + 'iml_2017/metadata_vectors/griml_2017_inventory_final_first_intervention.shp'
-    outtxt1 = workspace1 + 'iml_2017/metadata_vectors/generalstats.csv'
+    # #Define input file and location 
+    # workspace1 = '/home/pho/python_workspace/GrIML/other/'
+    # file1 = workspace1 + 'iml_2017/metadata_vectors/griml_2017_inventory_final_first_intervention.shp'
+    # outtxt1 = workspace1 + 'iml_2017/metadata_vectors/generalstats.csv'
+    
+    # general_stats(file1, outtxt1)
+    
     
-    general_stats(file1, outtxt1)
+    workspace1 = '/home/pho/python_workspace/GrIML/other/iml_2016-2023/final/checked/*IML-fv1.shp'
+    out_dir = '/home/pho/python_workspace/GrIML/other/iml_2016-2023/final/checked/stats/'
+    for f in list(glob.glob(workspace1)):
+        name = str(Path(f).stem)
+        year = str(Path(f).stem)[0:4]
+        out_file = str(Path(out_dir).joinpath(name+'_general_stats.txt'))
+        
+        print('\n\n'+str(Path(f).stem))
+        readout = general_stats(f, out_file)
+