From 9d988684cabc76c66074394b52f049e2236e5658 Mon Sep 17 00:00:00 2001 From: Cameron Scott Bodine Date: Mon, 24 Jun 2024 10:38:40 -0700 Subject: [PATCH] Test rectifying individual pings #120 --- src/class_portstarObj.py | 36 ++++-- src/class_rectObj.py | 55 +++++++-- src/class_sonObj.py | 30 ++++- src/main_readFiles.py | 259 ++++++++++++++++++++------------------- src/main_rectify.py | 22 +++- 5 files changed, 248 insertions(+), 154 deletions(-) diff --git a/src/class_portstarObj.py b/src/class_portstarObj.py index c871b48..6500158 100644 --- a/src/class_portstarObj.py +++ b/src/class_portstarObj.py @@ -348,7 +348,8 @@ def _createMosaicTransect(self, overview=True, threadCnt=cpu_count(), son=True, - maxChunk = 50): + maxChunk = 50, + cog=True): ''' Main function to mosaic exported rectified sonograms into a mosaic. If overview=True, overviews of the mosaic will be built, enhancing view @@ -384,6 +385,13 @@ def _createMosaicTransect(self, # maxChunk = 50 # Max chunks per mosaic. Limits each mosaic file size. self.imgsToMosaic = [] # List to store files to mosaic. + if cog: + chunkField = 'chunk_id' + else: + chunkField = 'chunk_id_2' + + print(chunkField) + if son: if self.port.rect_wcp: # Moscaic wcp sonograms if previousl exported self.port._loadSonMeta() @@ -393,7 +401,7 @@ def _createMosaicTransect(self, port = [] for name, group in df.groupby('transect'): - chunks = pd.unique(group['chunk_id']) + chunks = pd.unique(group[chunkField]) port_transect = [] for chunk in chunks: img_path = os.path.join(portPath, '*{}.tif'.format(chunk)) @@ -408,7 +416,7 @@ def _createMosaicTransect(self, star = [] for name, group in df.groupby('transect'): - chunks = pd.unique(group['chunk_id']) + chunks = pd.unique(group[chunkField]) star_transect = [] for chunk in chunks: img_path = os.path.join(starPath, '*{}.tif'.format(chunk)) @@ -427,12 +435,15 @@ def _createMosaicTransect(self, port = [] for name, group in df.groupby('transect'): - chunks = pd.unique(group['chunk_id']) + chunks = pd.unique(group[chunkField]) port_transect = [] for chunk in chunks: - img_path = os.path.join(portPath, '*{}.tif'.format(chunk)) - img = glob(img_path)[0] - port_transect.append(img) + try: + img_path = os.path.join(portPath, '*{}.tif'.format(chunk)) + img = glob(img_path)[0] + port_transect.append(img) + except: + pass port.append(port_transect) self.star._loadSonMeta() @@ -442,12 +453,15 @@ def _createMosaicTransect(self, star = [] for name, group in df.groupby('transect'): - chunks = pd.unique(group['chunk_id']) + chunks = pd.unique(group[chunkField]) star_transect = [] for chunk in chunks: - img_path = os.path.join(starPath, '*{}.tif'.format(chunk)) - img = glob(img_path)[0] - star_transect.append(img) + try: + img_path = os.path.join(starPath, '*{}.tif'.format(chunk)) + img = glob(img_path)[0] + star_transect.append(img) + except: + pass star.append(star_transect) srcToMosaic = [list(itertools.chain(*i)) for i in zip(port, star)] diff --git a/src/class_rectObj.py b/src/class_rectObj.py index a2c4482..e1e41c9 100644 --- a/src/class_rectObj.py +++ b/src/class_rectObj.py @@ -447,7 +447,8 @@ def _applyPosOffset(self, x_offset, y_offset): #=========================================== def _getRangeCoords(self, flip = False, - filt = 25): + filt = 25, + cog = True): ''' Humminbird SSS store one set geographic coordinates where each ping orriginates from (assuming GPS is located directly above sonar transducer). @@ -520,7 +521,7 @@ def _getRangeCoords(self, rotate *= -1 # Calculate ping bearing and normalize to range 0-360 - cog = True + # cog = False if cog: sDF[ping_bearing] = (sDF['cog']+rotate) % 360 else: @@ -586,7 +587,20 @@ def _getRangeCoords(self, ########################################## # Smooth and interpolate range coordinates - self._interpRangeCoords(filt) + if cog: + self._interpRangeCoords(filt) + else: + sDF = sDF[['record_num', 'chunk_id', 'ping_cnt', 'time_s', 'lons', 'lats', 'utm_es', 'utm_ns', 'instr_heading', 'cog', 'dep_m', 'range', 'range_lon', 'range_lat', 'range_e', 'range_n']].copy() + sDF.rename(columns={'lons': 'trk_lons', 'lats': 'trk_lats', 'utm_es': 'trk_utm_es', 'utm_ns': 'trk_utm_ns', 'cog': 'trk_cog', 'range_lat':'range_lats', 'range_lon':'range_lons', 'range_e':'range_es', 'range_n':'range_ns'}, inplace=True) + sDF['chunk_id_2'] = sDF.index.astype(int) + + ########################################### + # Overwrite Trackline_Smth_son.beamName.csv + # outCSV = os.path.join(self.metaDir, "Trackline_Smth_"+self.beamName+".csv") + outCSV = self.smthTrkFile + sDF.to_csv(outCSV, index=True, float_format='%.14f') + + # sys.exit() gc.collect() self._pickleSon() return #self @@ -1122,6 +1136,7 @@ def _exportCovShp(self, def _rectSonParallel(self, chunk, filt=50, + cog=True, wgs=False, son=True): ''' @@ -1209,7 +1224,14 @@ def _rectSonParallel(self, self._loadSonMeta() sonMetaAll = self.sonMetaDF - isChunk = sonMetaAll['chunk_id']==chunk + if cog: + isChunk = sonMetaAll['chunk_id']==chunk + else: + isChunk = sonMetaAll['chunk_id_2']==chunk + # next = sonMetaAll['chunk_id_2']==(chunk+1) + # isChunk = pd.concat([isChunk, next], ignore_index=True) + isChunk.iloc[chunk+1] = True + sonMeta = sonMetaAll[isChunk].reset_index() # Update class attributes based on current chunk @@ -1219,7 +1241,7 @@ def _rectSonParallel(self, if son: # Open image to rectify - self._getScanChunkSingle(chunk) + self._getScanChunkSingle(chunk, cog) else: # Rectifying substrate classification pass @@ -1234,6 +1256,9 @@ def _rectSonParallel(self, del self.shadowMask img = self.sonDat + # if not cog: + # # Zero out second ping + # img[:,1] = 0 # For each ping, we need the pixel coordinates where the sonar ## originates on the trackline, and where it terminates based on the @@ -1248,9 +1273,12 @@ def _rectSonParallel(self, # Create mask for filtering array. This makes fitting PiecewiseAffineTransform ## more efficient mask = np.zeros(len(pixAll), dtype=bool) # Create mask same size as pixAll - mask[0::filt] = 1 # Filter row coordinates - mask[1::filt] = 1 # Filter column coordinates - mask[-2], mask[-1] = 1, 1 # Make sure we keep last row/col coordinates + if cog: + mask[0::filt] = 1 # Filter row coordinates + mask[1::filt] = 1 # Filter column coordinates + mask[-2], mask[-1] = 1, 1 # Make sure we keep last row/col coordinates + else: + mask[:] = 1 # Filter pix pix = pixAll[mask] @@ -1262,7 +1290,16 @@ def _rectSonParallel(self, # Open smoothed trackline/range extent file trkMeta = pd.read_csv(trkMetaFile) - trkMeta = trkMeta[trkMeta['chunk_id']==chunk].reset_index(drop=False) # Filter df by chunk_id + if cog: + trkMeta = trkMeta[trkMeta['chunk_id']==chunk].reset_index(drop=False) # Filter df by chunk_id + else: + # trkMeta = trkMeta[trkMeta['chunk_id_2']==chunk].reset_index(drop=False) + # next = trkMeta[trkMeta['chunk_id_2']==chunk+1].reset_index(drop=False) + # trkMeta = pd.concat([trkMeta, next], ignore_index=True) + isChunk = trkMeta['chunk_id_2']==chunk + isChunk.iloc[chunk+1] = True + trkMeta = trkMeta[isChunk].reset_index(drop=False) + pix_m = self.pixM # Get pixel size # Get range (outer extent) coordinates [xR, yR] to transposed numpy arrays diff --git a/src/class_sonObj.py b/src/class_sonObj.py index 695638f..55e92c7 100644 --- a/src/class_sonObj.py +++ b/src/class_sonObj.py @@ -1909,6 +1909,7 @@ def _doSpdCor(self, chunk, lbl_set=1, spdCor=1, maxCrop=0, son=True, integer=Tru # ====================================================================== def _getScanChunkSingle(self, chunk, + cog=True, filterIntensity = False, remWater = False): ''' @@ -1946,7 +1947,11 @@ def _getScanChunkSingle(self, sonMetaAll = pd.read_csv(self.sonMetaFile) # Filter df by chunk - isChunk = sonMetaAll['chunk_id']==chunk + if cog: + isChunk = sonMetaAll['chunk_id']==chunk + else: + isChunk = sonMetaAll['chunk_id_2']==chunk + isChunk.iloc[chunk+1] = True sonMeta = sonMetaAll[isChunk].reset_index() # Update class attributes based on current chunk @@ -1994,6 +1999,29 @@ def _getChunkID(self): del self.sonMetaDF, df return chunks + + # ====================================================================== + def _getChunkID_Update(self): + ''' + Utility to load unique chunk ID's from son obj and return in a list + ''' + + # Load son metadata csv to df + self._loadSonMeta() + + # # Get unique chunk id's + # df = self.sonMetaDF.groupby(['chunk_id', 'index']).size().reset_index().rename(columns={0:'count'}) + # chunks = pd.unique(df['chunk_id']).astype(int) + + # Use index as chunk id + df = self.sonMetaDF + chunks = df.index.values.astype(int) + + df['chunk_id_2'] = chunks + self._saveSonMetaCSV(df) + + del self.sonMetaDF, df + return chunks # ====================================================================== def _addZero(self, chunk): diff --git a/src/main_readFiles.py b/src/main_readFiles.py index 1fe2112..473a899 100644 --- a/src/main_readFiles.py +++ b/src/main_readFiles.py @@ -272,6 +272,11 @@ def read_master_func(logfilename='', ## and access sonar data. We will use the first sonar beam to make an ## initial sonar object, then create a copy for each beam. tempC = float(tempC)/10 # Divide temperature by 10 + + # cog = False + # if not cog: + # nchunk= 1 + son = sonObj(sonFiles[0], humFile, projDir, tempC, nchunk) # Initialize sonObj instance from first sonar beam ####################################### @@ -507,7 +512,7 @@ def read_master_func(logfilename='', # Get metadata for each beam in parallel. if len(toProcess) > 0: - r = Parallel(n_jobs= np.min([len(toProcess), threadCnt]), verbose=10)(delayed(son._getSonMeta)() for son in toProcess) + # r = Parallel(n_jobs= np.min([len(toProcess), threadCnt]), verbose=10)(delayed(son._getSonMeta)() for son in toProcess) # Store pix_m in object for son, pix_m in zip(sonObjs, r): son.pixM = pix_m # Sonar instrument pixel resolution @@ -898,137 +903,137 @@ def read_master_func(logfilename='', # Use AOI to clip metadata to facilitate processing multiple transects in a ## single sonar recording. - if aoi: - - print('\n\nProcessing AOI...') - - # # # For extracting value from nested dictionaries and lists - # def getPolyCoords(nested_dict, value): - # # polys = [] - # for k, v in nested_dict.items(): - # if k == value: - # return v - # elif hasattr(v, 'items'): - # p = getPolyCoords(v, value) - # if p is not None and p: - # # return p - # polys.append(p) - # elif isinstance(v, list): - # for i in v: - # if hasattr(i, 'items'): - # p = getPolyCoords(i, value) - # if p is not None and p: - # # return p - # polys.append(p) - # return polys - - # If .plan file (from Hydronaulix) - if os.path.basename(aoi.split('.')[-1]) == 'plan': - with open(aoi, 'r', encoding='utf-8') as f: - f = json.load(f) - # Find 'polygon' coords in nested json - # polys = [] - # poly_coords = getPolyCoords(f, 'polygon') - # print(poly_coords) - - f = f['mission'] - f = f['items'] - poly_coords = [] - for i in f: - for k, v in i.items(): - if k == 'polygon': - poly_coords.append(v) - - aoi_poly_all = gpd.GeoDataFrame() - - for poly in poly_coords: + # if aoi: + + # print('\n\nProcessing AOI...') + + # # # # For extracting value from nested dictionaries and lists + # # def getPolyCoords(nested_dict, value): + # # # polys = [] + # # for k, v in nested_dict.items(): + # # if k == value: + # # return v + # # elif hasattr(v, 'items'): + # # p = getPolyCoords(v, value) + # # if p is not None and p: + # # # return p + # # polys.append(p) + # # elif isinstance(v, list): + # # for i in v: + # # if hasattr(i, 'items'): + # # p = getPolyCoords(i, value) + # # if p is not None and p: + # # # return p + # # polys.append(p) + # # return polys + + # # If .plan file (from Hydronaulix) + # if os.path.basename(aoi.split('.')[-1]) == 'plan': + # with open(aoi, 'r', encoding='utf-8') as f: + # f = json.load(f) + # # Find 'polygon' coords in nested json + # # polys = [] + # # poly_coords = getPolyCoords(f, 'polygon') + # # print(poly_coords) + + # f = f['mission'] + # f = f['items'] + # poly_coords = [] + # for i in f: + # for k, v in i.items(): + # if k == 'polygon': + # poly_coords.append(v) + + # aoi_poly_all = gpd.GeoDataFrame() + + # for poly in poly_coords: - # Extract coordinates - lat_coords = [i[0] for i in poly] - lon_coords = [i[1] for i in poly] - - polygon_geom = Polygon(zip(lon_coords, lat_coords)) - aoi_poly = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polygon_geom]) - - aoi_poly_all = pd.concat([aoi_poly_all, aoi_poly], ignore_index=True) - - # If shapefile - elif os.path.basename(aoi.split('.')[-1]) == 'shp': - aoi_poly_all = gpd.read_file(aoi) - - else: - print(os.path.basename, ' is not a valid aoi file type.') - sys.exit() - - # Reproject to utm - epsg = int(sonObjs[0].humDat['epsg'].split(':')[-1]) - aoi_poly = aoi_poly_all.to_crs(crs=epsg) - aoi_poly = aoi_poly.dissolve() - - # Buffer aoi - if os.path.basename(aoi.split('.')[-1]) == 'plan': - buf_dist = 0.5 - aoi_poly['geometry'] = aoi_poly.geometry.buffer(buf_dist) - - # Save aoi - aoi_dir = os.path.join(sonObjs[0].projDir, 'aoi') - aoiOut = os.path.basename(sonObjs[0].projDir) + '_aoi.shp' - if not os.path.exists(aoi_dir): - os.makedirs(aoi_dir) - - aoiOut = os.path.join(aoi_dir, aoiOut) - aoi_poly.to_file(aoiOut) - - # Iterate each son file, clip with aoi, and save - for son in sonObjs: - son._loadSonMeta() - sonDF = son.sonMetaDF - - # Convert to geodataframe - epsg = int(son.humDat['epsg'].split(':')[-1]) - sonDF = gpd.GeoDataFrame(sonDF, geometry=gpd.points_from_xy(sonDF.e, sonDF.n), crs=epsg) - - # Clip with aoi - sonDF = gpd.clip(sonDF, aoi_poly, keep_geom_type=False) - sonDF = sonDF.sort_index() - - # Make transects from consective pings using dataframe index - idx = sonDF.index.values - transect_groups = np.split(idx, np.where(np.diff(idx) != 1)[0]+1) - - # Assign transect - transect = 0 - for t in transect_groups: - sonDF.loc[sonDF.index>=t[0], 'transect'] = transect - transect += 1 - - # Set chunks - lastChunk = 0 - newChunk = [] - for name, group in sonDF.groupby('transect'): - - if (len(group)%nchunk) != 0: - rdr = nchunk-(len(group)%nchunk) - chunkCnt = int(len(group)/nchunk) - chunkCnt += 1 - else: - rdr = False - chunkCnt = int(len(group)/nchunk) - - chunks = np.arange(chunkCnt) + lastChunk - chunks = np.repeat(chunks, nchunk) + # # Extract coordinates + # lat_coords = [i[0] for i in poly] + # lon_coords = [i[1] for i in poly] + + # polygon_geom = Polygon(zip(lon_coords, lat_coords)) + # aoi_poly = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polygon_geom]) + + # aoi_poly_all = pd.concat([aoi_poly_all, aoi_poly], ignore_index=True) + + # # If shapefile + # elif os.path.basename(aoi.split('.')[-1]) == 'shp': + # aoi_poly_all = gpd.read_file(aoi) + + # else: + # print(os.path.basename, ' is not a valid aoi file type.') + # sys.exit() + + # # Reproject to utm + # epsg = int(sonObjs[0].humDat['epsg'].split(':')[-1]) + # aoi_poly = aoi_poly_all.to_crs(crs=epsg) + # aoi_poly = aoi_poly.dissolve() + + # # Buffer aoi + # if os.path.basename(aoi.split('.')[-1]) == 'plan': + # buf_dist = 0.5 + # aoi_poly['geometry'] = aoi_poly.geometry.buffer(buf_dist) + + # # Save aoi + # aoi_dir = os.path.join(sonObjs[0].projDir, 'aoi') + # aoiOut = os.path.basename(sonObjs[0].projDir) + '_aoi.shp' + # if not os.path.exists(aoi_dir): + # os.makedirs(aoi_dir) + + # aoiOut = os.path.join(aoi_dir, aoiOut) + # aoi_poly.to_file(aoiOut) + + # # Iterate each son file, clip with aoi, and save + # for son in sonObjs: + # son._loadSonMeta() + # sonDF = son.sonMetaDF + + # # Convert to geodataframe + # epsg = int(son.humDat['epsg'].split(':')[-1]) + # sonDF = gpd.GeoDataFrame(sonDF, geometry=gpd.points_from_xy(sonDF.e, sonDF.n), crs=epsg) + + # # Clip with aoi + # sonDF = gpd.clip(sonDF, aoi_poly, keep_geom_type=False) + # sonDF = sonDF.sort_index() + + # # Make transects from consective pings using dataframe index + # idx = sonDF.index.values + # transect_groups = np.split(idx, np.where(np.diff(idx) != 1)[0]+1) + + # # Assign transect + # transect = 0 + # for t in transect_groups: + # sonDF.loc[sonDF.index>=t[0], 'transect'] = transect + # transect += 1 + + # # Set chunks + # lastChunk = 0 + # newChunk = [] + # for name, group in sonDF.groupby('transect'): + + # if (len(group)%nchunk) != 0: + # rdr = nchunk-(len(group)%nchunk) + # chunkCnt = int(len(group)/nchunk) + # chunkCnt += 1 + # else: + # rdr = False + # chunkCnt = int(len(group)/nchunk) + + # chunks = np.arange(chunkCnt) + lastChunk + # chunks = np.repeat(chunks, nchunk) - if rdr: - chunks = chunks[:-rdr] + # if rdr: + # chunks = chunks[:-rdr] - newChunk += list(chunks) - lastChunk = chunks[-1] + 1 - del chunkCnt + # newChunk += list(chunks) + # lastChunk = chunks[-1] + 1 + # del chunkCnt - sonDF['chunk_id'] = newChunk + # sonDF['chunk_id'] = newChunk - son._saveSonMetaCSV(sonDF) - son._cleanup() + # son._saveSonMetaCSV(sonDF) + # son._cleanup() ############################################################################ diff --git a/src/main_rectify.py b/src/main_rectify.py index 49c0798..b186f75 100644 --- a/src/main_rectify.py +++ b/src/main_rectify.py @@ -70,6 +70,7 @@ def rectify_master_func(logfilename='', smthDep=0, adjDep=0, pltBedPick=False, + cog=False, rect_wcp=False, rect_wcr=False, son_colorMap='Greys', @@ -354,9 +355,14 @@ def rectify_master_func(logfilename='', ############################################################################ # Calculate range extent coordinates # ############################################################################ + # cog=True + start_time = time.time() - print("\nCalculating, smoothing, and interpolating range extent coordinates...") - Parallel(n_jobs= np.min([len(portstar), threadCnt]), verbose=10)(delayed(son._getRangeCoords)(flip, filterRange) for son in portstar) + if cog: + print("\nCalculating, smoothing, and interpolating range extent coordinates...") + else: + print("\nCalculating range extent coordinates from vessel heading...") + Parallel(n_jobs= np.min([len(portstar), threadCnt]), verbose=10)(delayed(son._getRangeCoords)(flip, filterRange, cog) for son in portstar) print("Done!") print("Time (s):", round(time.time() - start_time, ndigits=1)) gc.collect() @@ -406,7 +412,11 @@ def rectify_master_func(logfilename='', son.outDir = os.path.join(son.projDir, son.beamName) # Get chunk id's - chunks = son._getChunkID() + if cog: + chunks = son._getChunkID() + else: + chunks = son._getChunkID_Update() + chunks = chunks[:-1] # Load sonMetaDF son._loadSonMeta() @@ -416,9 +426,9 @@ def rectify_master_func(logfilename='', print('\n\tExporting', len(chunks), 'GeoTiffs for', son.beamName) # for i in chunks: - # son._rectSonParallel(i, filter, wgs=False) + # son._rectSonParallel(i, filter, cog, wgs=False) # sys.exit() - Parallel(n_jobs= np.min([len(chunks), threadCnt]), verbose=10)(delayed(son._rectSonParallel)(i, filter, wgs=False) for i in chunks) + Parallel(n_jobs= np.min([len(chunks), threadCnt]), verbose=10)(delayed(son._rectSonParallel)(i, filter, cog, wgs=False) for i in chunks) son._cleanup() gc.collect() printUsage() @@ -457,7 +467,7 @@ def rectify_master_func(logfilename='', if not aoi: psObj._createMosaic(mosaic, overview, threadCnt, son=True, maxChunk=mosaic_nchunk) else: - psObj._createMosaicTransect(mosaic, overview, threadCnt, son=True, maxChunk=mosaic_nchunk) + psObj._createMosaicTransect(mosaic, overview, threadCnt, son=True, maxChunk=mosaic_nchunk, cog=cog) print("Done!") print("Time (s):", round(time.time() - start_time, ndigits=1)) del psObj