Merging functions from seabed_tools to geography

- Deleting old functions in seabed_tools that are unused nowadays (getPlotBounds, getCoast, getDepthFromBathymetryMesh)
- - convertLatLong2Meters - old function that has now been updated in geography
- - processASC - basically the old version of getMapBathymetry - also in geography
- Merging writeBathymetryFile functions between seabed_tools and geography
- processGeoTiff now moved to geography since it deals with processing geographic-related information

Author: shousner

famodel/geography.py

@@ -337,14 +337,21 @@ def writeBathymetryFile(moorpy_bathymetry_filename, bathXs, bathYs, bath_depths,
     '''Write a MoorDyn/MoorPy-style bathymetry text file based on provided
     x and y grid line values and a 2D array of depth values.'''
 
+    # open the file
     f = open(os.path.join(os.getcwd(), moorpy_bathymetry_filename), 'w')
+
     f.write('--- MoorPy Bathymetry Input File ---\n')
     f.write(f'nGridX {len(bathXs)}\n')
     f.write(f'nGridY {len(bathYs)}\n')
     f.write(f'      ')
+
+    # x-coordinates
     for ix in range(len(bathXs)):
         f.write(f'{bathXs[ix]:.2f} ')
     f.write('\n')
+    #f.write(" ".join(map(str, grid_x)) + "\n")  # old version from seabed_tools
+
+    # y-coordintes
     for iy in range(len(bathYs)):
         f.write(f'{bathYs[iy]:.2f} ')
         for id in range(len(bath_depths[iy])):
@@ -353,6 +360,11 @@ def writeBathymetryFile(moorpy_bathymetry_filename, bathXs, bathYs, bath_depths,
             else:
                 f.write(f'{bath_depths[iy,id]:8.3f} ')
         f.write('\n')
+    #for i, y in enumerate(grid_y):         # alternative writing version
+        #row = [y] + list(grid_depth[i, :])
+        #f.write(" ".join(map(str, row)) + "\n")
+    
+    # close the file
     f.close()
 
 
@@ -407,6 +419,94 @@ def getLeaseAndBathymetryInfo(lease_name, bathymetry_file, bath_ncols=100, bath_
 
 
 
+
+def processGeotiff(filename, lat, lon, outfilename="processGeotiff.txt", **kwargs):
+    '''Process a geotiff file containing bathymetry (or other info)
+    and convert into a rectangular bathymetry grid in units of m relative to 
+    the project reference point.
+    
+    Parameters
+    ----------
+    filename : string
+        Path and name of geotiff file to load.
+    lat : float
+        lattitude of reference point to use for array y grid
+    long : float
+        lattitude of reference point to use for array x grid
+    outfilename : string, optional
+        If provided, writes a MoorDyn/MoorPy style bathymetry file
+    kwargs : dict
+        Optional extra arguments that will be relayed to convertBathymetry2Meters,
+        can be used to specify desired x and y grid coordinates.
+
+    Returns
+    -------
+    Xs : array
+        x values of grid points [m]
+    Ys : array
+        y values of grid points [m]
+    depths  : 2D array
+        water depth grid (positive down) [m]
+    '''
+
+    import rasterio
+    import rasterio.plot
+
+    tiff = rasterio.open(filename)  # load the geotiff file
+    
+    #rasterio.plot.show(tiff)  # plot it to see that it works
+    
+    # note: a CRS is stored with the geotiff, accessible with tiff.crs
+    
+    # Get lattitude and longitude grid values
+    #_, longs = rasterio.transform.xy(tiff.transform, range(tiff.height),0)
+    #lats, _  = rasterio.transform.xy(tiff.transform, 0, range(tiff.width-1,-1,-1))
+    height, width = tiff.shape
+    cols, rows = np.meshgrid(np.arange(width), np.arange(height))
+    longs_mesh, lats_mesh = rasterio.transform.xy(tiff.transform, rows, cols)
+    longs = np.array(longs_mesh)[0,:]
+    lats = np.flip(np.array(lats_mesh)[:,0])
+    # lats data provided from top left corner, i.e., latitudes are descending. It seems that the following interpolation functions (getDepthFromBathymetry)
+    # can only work if latitudes start small and increase, meaning that the first latitude entry has to be the bottom left corner
+    
+    # Depth values in numpy array form
+    depths = -tiff.read(1)
+    depths = np.flipud(depths)
+    # because the interpolation functions require the latitude array to go from small to big (i.e., the bottom left corner), we need to flip the depth matrix to align
+    # it will all get sorted out later to what it should be geographically when plotting in MoorPy
+    
+    
+    # extract the coordinate reference systems needed (pyproj CRS objects)
+    latlong_crs = tiff.crs
+    #target_crs = getTargetCRS(lon, lat)     # should be UTM 10N for Humboldt/California coast
+    target_crs = getCustomCRS(lon, lat)     # get custom CRS centered around the lat/long point you want
+    
+    # get the centroid/reference location in lat/long coordinates
+    centroid = (lon, lat)
+
+    # get the centroid/reference location in target_crs coordinates
+    #centroid_utm = (lon, lat)
+    _, _, centroid_utm = convertLatLong2Meters([centroid[0]], [centroid[1]], centroid, latlong_crs, target_crs, return_centroid=True)
+    
+    # set the number of rows and columns to use in the MoorPy bathymetry file
+    ncols = 100
+    nrows = 100
+    
+    # convert bathymetry to meters
+    bath_xs, bath_ys, bath_depths = convertBathymetry2Meters(longs, lats, depths, 
+                                                             centroid, centroid_utm, 
+                                                             latlong_crs, target_crs, 
+                                                             ncols, nrows, **kwargs)
+    # export to MoorPy-readable file
+    writeBathymetryFile(outfilename, bath_xs, bath_ys, bath_depths)
+
+    return bath_xs, bath_ys, bath_depths
+
+
+
+
+
+
 def getSoilType(x, y, centroid, latlong_crs, custom_crs, soil_file):
     """Function to return the name of the soil below a specific x/y coordinate by creating shapely polygons based on the shapefile data.
     It loops through all polygons in the shapefile and if the x/y position is contained in that polygon, it returns the soil of that polygon."""