Towards cleaning up the geography module and seabed_tools

shousner · shousner · commit 0f2bc42abbd0 · 2025-10-30T22:37:49.000-06:00
- Deleted projectAlongSeabed function since that has now moved to main project.py file
- Deleted all of the old commented out code (that was used for plots in the 2023 shared anchor LCOE report)
- Deleted the plot3d function because this is also outdated from the new project.plot3d
- Still some more ideas to be addressed in merging in similar ideas from seabed_tools
diff --git a/famodel/geography.py b/famodel/geography.py
@@ -520,377 +520,7 @@ def getSoilGrid(centroid, latlong_crs, custom_crs, soil_file, nrows=100, ncols=1
 
 
         
-def plot3d(lease_xs, lease_ys, bathymetryfilename, area_on_bath=False, args_bath={}, xbounds=None, ybounds=None, zbounds=None):
-    '''Plot aspects of the Project object in matplotlib in 3D'''
-
-    # organize the bathymetry arguments
-    if len(args_bath)==0:
-        args_bath = {'zlim':[-3200,500], 'cmap':'gist_earth'}
- 
-    fig = plt.figure(figsize=(6,4))
-    ax = plt.axes(projection='3d')
-
-    if xbounds != None:
-        ax.set_xlim(xbounds[0], xbounds[1])
-    if ybounds != None:
-        ax.set_ylim(ybounds[0], ybounds[1])
-    if zbounds != None:
-        ax.set_zlim(zbounds[0], zbounds[1])
-
-    # plot the lease area in a red color, if desired
-    ax.plot(lease_xs, lease_ys, np.zeros(len(lease_xs)), color='r', zorder=100)
-    
-    # plot the bathymetry in matplotlib using a plot_surface
-
-    # !!!! include option to plot entire bathymetry file or not
-
-    if isinstance(bathymetryfilename, str):
-        bathGrid_Xs, bathGrid_Ys, bathGrid = sbt.readBathymetryFile(bathymetryfilename)         # parse through the MoorDyn/MoorPy-formatted bathymetry file
-        X, Y = np.meshgrid(bathGrid_Xs, bathGrid_Ys)                                    # create a 2D mesh of the x and y values
-        bath = ax.plot_surface(X, Y, -bathGrid, rstride=1, cstride=1,
-                               vmin=args_bath['zlim'][0], vmax=args_bath['zlim'][1], 
-                               cmap=args_bath['cmap'])
-    
-    '''
-    # plot the project boundary
-    if boundary:
-        ax.plot(self.boundaryXs, self.boundaryYs, np.zeros(len(self.boundaryXs)), color='b', zorder=100, alpha=0.5)
-    '''
-        
-    # plot the projection of the lease area bounds on the seabed, if desired
-    if area_on_bath:
-        lease_zs = projectAlongSeabed(lease_xs, lease_ys, bathGrid_Xs, bathGrid_Ys, bathGrid)
-        ax.plot(lease_xs, lease_ys, -lease_zs, color='tab:blue', zorder=10, alpha=0.5)
-
-
-    set_axes_equal(ax)
-    ax.axis('off')
-
-    return fig, ax
-
-
-
-def projectAlongSeabed(x, y, bathXs, bathYs, bath_depths):
-    '''Project a set of x-y coordinates along a seabed surface (grid),
-    returning the corresponding z coordinates.'''
-    
-    if len(x) == len(y):
-        n = len(x)        
-        z = np.zeros(n)   # z coordinate of each point [m]
-        a = np.zeros(n)   # could also do slope (dz/dh)
-        for i in range(n):
-            z[i], nvec = sbt.getDepthFromBathymetry(x[i], y[i], bathXs, bathYs, bath_depths)
-    
-    else:
-        z = np.zeros([len(y), len(x)])
-        for i in range(len(y)):
-            for j in range(len(x)):
-                z[i,j], nvec = sbt.getDepthFromBathymetry(x[j], y[i], bathXs_mesh, bathYs_mesh, bath_depths)
-        
-    return z
-
-
-
-
-
-
-
-
-        
-
-
-    """
-        if self.lat0 != 0 and self.lon0 != 0:
-                # set the centroid of the project and save in a GeoDataFrame
-                self.centroid = (self.lon0, self.lat0)
-                self.gdf = gpd.GeoDataFrame({'type':'centroid', 'geometry': [Point(self.centroid)]}, crs=self.latlong_crs)
-                # set the target coordinate reference system (CRS) that will switch between regular lat/long system and "meters from centroid" system, based on centroid location
-                information
-                            # extract the numeric code
-                self.  # save the target CRS (UTM 10N = 32610)
-        
-
-        gdf_leases = gpd.GeoDataFrame({'type': 'lease_area', 'geometry': lease_area.geometry}, crs=getLatLongCRS() )
-
-        gdf = add2geodataframe(gdf0, gdf_leases)
-
-    # convert lat/long boundary/lease area points to meters away from the centroid
-            self.lease_xs, self.lease_ys = self.convertLatLong2Meters(self.area_longs, self.area_lats, self.centroid)
-
-
-            if which_centroid=='lease_area':
-                # make a blank copy of the gdf to switch to the target CRS to get the accurate centroid
-                gdf_utm = self.gdf.copy().to_crs(self.target_crs)
-                centroid_utm = (gdf_utm.geometry.centroid.values.x[0], gdf_utm.geometry.centroid.values.y[0])
-                gdf_centroid = gpd.GeoDataFrame({'type':'centroid', 'geometry': [Point(centroid_utm)]}, crs=self.target_crs)
-            else:
-                # assuming the input centroid is in a long/lat pair
-                gdf_centroid = gpd.GeoDataFrame({'type':'centroid', 'geometry': [Point(which_centroid)]}, crs=self.latlong_crs)
-                gdf_centroid.to_crs(self.target_crs)
-            
-            self.centroid_utm = (gdf_centroid.geometry.values.x[0], gdf_centroid.geometry.values.y[0])
-            self.centroid = (gdf_centroid.to_crs(self.latlong_crs).geometry.values.x[0], gdf_centroid.to_crs(self.latlong_crs).geometry.values.y[0]) # assume centroid is focal point of Project
-
-            def initialize_geodataframe(crs, columns=['type','geometry']):
-
-        gdf = gpd.GeoDataFrame(columns=columns, geometry='initialization', crs=crs)
-
-        return gdf
-
-    def add2geodataframe(gdf_to_add_to, gdf_to_add):
-        
-        # check to make sure they have the same columns and CRS
-
-        gdf_new = pd.concat([gdf_to_add_to, gdf_to_add])
-
-        return gdf_new
-
-
-
-
-
-        # TODO
-
-        # reference entire CA bathymetry file (maybe)
-        # plot2d method and a plotGDF method (with bathymetry in the geodataframe using contours)
-
-        # add the coastline
-        #xcoast, ycoast = sbt.getCoast(self.Xs, self.Ys, self.depths)
-        #ax.plot(xcoast, ycoast, np.zeros(len(self.Ys)), color='k', zorder=100)
-
-        # need to fix up bounds
-        #xbmin, xbmax = sbt.getPlotBounds(self.longs_bath, self.centroid, long=True)
-        #ybmin, ybmax = sbt.getPlotBounds(self.lats_bath, self.centroid, long=False)
-
-        plt.show()
-
-
-
-
-
-
-
-
-    # METHODS USED SPECIFICALLY FOR DEEPFARM LCOE ANALYSIS
-    
-    def addMap2GDF(self, filename='', states=None):
-        '''function to include a shapefile of a provided map'''
-
-        # read in the provided filename to add to the geodataframe
-        usa = gpd.read_file(filename)
-        # internal list of states, in order, to go with the default U.S. states shapefile
-        statenamelist = ['Maryland','Iowa','Delaware','Ohio','Pennsylvania','Nebraska','Washington','Puerto Rico','Alabama','Arkansas','New Mexico',            # 0-10
-                         'Texas','California','Kentucky','Georgia','Wisconsin','Oregon','Missouri','Virginia','Tennessee','Louisiana','New York',               # 11-21
-                         'Michigan','Idaho','Florida','Alaska','Illinois','Montana','Minnesota','Indiana','Massachusetts','Kansas','Nevada','Vermont',          # 22-33
-                         'Connecticut','New Jersey','Washington D.C.','North Carolina','Utah','North Dakota','South Carolina','Mississippi','Colorado',         # 34-42
-                         'South Dakota','Oklahoma','Wyoming','West Virginia','Maine','Hawaii','New Hampshire','Arizona','Rhode Island']                         # 43-51
-        # insert names of the states into the new gdf
-        usa.insert(0, 'type', statenamelist)
-        # set the CRS of the USA pdf to the right CRS
-        usa.set_crs(crs="EPSG:4326", inplace=True)
-        self.usa = usa
-        
-        for state in states:
-            state_gs = usa.loc[usa['type']==state]
-            self.gdf = pd.concat([self.gdf, state_gs])
-    
-
-
-    def setFarmLayout(self, style='grid', nrows=10, ncols=10, turbine_spacing=2000, nOSS=2):
-        
-        if style=='grid':
-            # for now, this is very custom code specific to the DeepFarm project
-            farmxspacing = (nrows-1)*turbine_spacing
-            farmyspacing = (ncols-1)*turbine_spacing
-
-            turbine_distances_from_centroid = []
-            oss_distances_from_centroid = []
-            for j in reversed(range(ncols)):
-                for i in range(nrows):
-                    xpos = -(farmxspacing/2)+(i*turbine_spacing)
-                    ypos = -(farmyspacing/2)+(j*turbine_spacing)
-                    turbine_distances_from_centroid.append((xpos, ypos))
-            
-            # add positions of two offshore substations (OSSs)
-            oss_distances_from_centroid.append((11000.0, 5000.0))
-            oss_distances_from_centroid.append((11000.0, -5000.0))
-        
-        if style=='shared':
-            turbine_xspacing = np.sqrt(2000**2-1000**2)
-            turbine_yspacing = 2000
-            farmxspacing = turbine_xspacing*9
-            farmyspacing = turbine_yspacing*9
-
-            turbine_distances_from_centroid = []
-            oss_distances_from_centroid = []
-            for j in reversed(range(ncols)):
-                for i in range(nrows):
-                    xpos = -(farmxspacing/2)+(i*turbine_xspacing)
-                    ypos = -(farmyspacing/2)+(j*turbine_yspacing) - 1000*np.sin(np.radians(30)) + 1000*(i%2)
-                    turbine_distances_from_centroid.append((xpos, ypos))
-            
-            # add positions of two offshore substations (OSSs)
-            oss_distances_from_centroid.append((5.5*turbine_xspacing, 2.0*turbine_yspacing+1000*np.sin(np.radians(30))))
-            oss_distances_from_centroid.append((5.5*turbine_xspacing, -2.5*turbine_yspacing-1000*np.sin(np.radians(30))))
-        
-        if style=='small-shared':
-            turbine_xspacing = np.sqrt(2000**2-1000**2)
-            turbine_yspacing = 2000
-            farmxspacing = turbine_xspacing*1
-            farmyspacing = turbine_yspacing*2
-
-            turbine_distances_from_centroid = []
-            oss_distances_from_centroid = []
-            for j in reversed(range(3)):
-                for i in range(2):
-                    xpos = -(farmxspacing/2)+(i*turbine_xspacing)
-                    ypos = -(farmyspacing/2)+(j*turbine_yspacing) - 1000*np.sin(np.radians(30)) + 1000*(i%2)
-                    turbine_distances_from_centroid.append((xpos, ypos))
-            
-            # add positions of two offshore substations (OSSs)
-            oss_distances_from_centroid.append((-0.5*turbine_xspacing, 2.0*turbine_yspacing-1000*np.sin(np.radians(30))))
-        
-
-        # create a copy of the global gdf and transform it into the easting/northing coordinate reference system
-        gdf_utm = self.gdf.copy().to_crs(self.target_crs)
-        xcentroid = gdf_utm.loc[gdf_utm['type']=='centroid'].centroid.x[0]
-        ycentroid = gdf_utm.loc[gdf_utm['type']=='centroid'].centroid.y[0]
-
-        # create shapely Point objects of the turbine positions relative to the centroid, in the UTM CRS
-        turbine_geoms = []
-        for i,(x,y) in enumerate(turbine_distances_from_centroid):
-            turbine_geoms.append( Point(xcentroid + x, ycentroid + y) )
-        
-        oss_geoms = []
-        for i,(x,y) in enumerate(oss_distances_from_centroid):
-            oss_geoms.append( Point(xcentroid + x, ycentroid + y) )
-        
-        # make a new gdf to put the turbine data together
-        turbine_gdf = gpd.GeoDataFrame({'type': ['turbine']*len(turbine_geoms), 'geometry': turbine_geoms}, crs=self.target_crs)
-        # make a new gdf to put the substation data together
-        oss_gdf = gpd.GeoDataFrame({'type': 'substation', 'geometry': oss_geoms}, crs=self.target_crs)
-        # merge these two geodataframes together into one (best way I can find to "add" rows to a dataframe; ignoring index makes all indices a different number)
-        turbine_gdf = pd.concat([turbine_gdf, oss_gdf], ignore_index=True)
-
-        # convert the turbine/oss coordinates back to regular latitude/longitude (EPSG: 4326)
-        turbine_gdf = turbine_gdf.to_crs('EPSG:4326')
-
-        # add the turbine gdf to the global gdf
-        self.gdf = pd.concat([self.gdf, turbine_gdf], ignore_index=True)
-
-        # add local variables in this method to the turbine_gdf to be used later (but don't need for the global gdf)
-        turbine_gdf.insert(2, 'easting_northing_geometry', turbine_geoms + oss_geoms)
-        turbine_gdf.insert(3, 'meters_from_centroid', turbine_distances_from_centroid + oss_distances_from_centroid)
-
-        # save this new turbine_gdf for future use
-        self.turbine_gdf = turbine_gdf
-
-        # make a layout CSV (used for WHaLE/WAVES)
-        self.makeLayoutCSV()
-
-    
-
-    def makeLayoutCSV(self, filename='layout_test.csv'):
-
-        turbine_longs = [point.coords[0][0] for point in self.turbine_gdf.geometry]
-        turbine_lats = [point.coords[0][1] for point in self.turbine_gdf.geometry]
-
-        self.turbine_gdf.insert(2, 'longitude', turbine_longs)
-        self.turbine_gdf.insert(3, 'latitude', turbine_lats)
-
-        turbine_eastings = [point.coords[0][0] for point in self.turbine_gdf['easting_northing_geometry']]
-        turbine_northings = [point.coords[0][1] for point in self.turbine_gdf['easting_northing_geometry']]
-
-        #self.turbine_gdf.insert(5, 'easting', turbine_eastings)
-        #self.turbine_gdf.insert(6, 'northing', turbine_northings)
-
-        turbine_x_from_centroid = [point[0] for point in self.turbine_gdf['meters_from_centroid']]
-        turbine_y_from_centroid = [point[1] for point in self.turbine_gdf['meters_from_centroid']]
-
-        self.turbine_gdf.insert(5, 'easting', turbine_x_from_centroid)
-        self.turbine_gdf.insert(6, 'northing', turbine_y_from_centroid)
-
-        self.turbine_gdf.insert(8, 'floris_x', turbine_x_from_centroid)
-        self.turbine_gdf.insert(9, 'floris_y', turbine_y_from_centroid)
-
-        columns = ['type', 'longitude', 'latitude', 'easting', 'northing', 'floris_x', 'floris_y']
-        df = pd.DataFrame(self.turbine_gdf)
-        df.to_csv(filename, columns=columns)
-
-
-    def plotGDF(self, kwargs):
-        '''2D map-like plot'''
-        
-        if 'centroid' in kwargs:
-            centroid_settings = kwargs['centroid']
-            if 'label' in centroid_settings:
-                centroid_label = 'centroid'
-        if 'map' in kwargs:
-            map_settings = kwargs['map']
-        if 'farm' in kwargs:
-            farm_settings = kwargs['farm']
-
-        fig, ax = plt.subplots(1,1)
-
-        if 'centroid' in kwargs:
-            self.gdf.loc[self.gdf['type']=='centroid'].plot(ax=ax, color=centroid_settings['color'], label=centroid_label)
-
-        if 'boundary' in kwargs:
-            map_boundary = self.gdf.loc[self.gdf['type']=='California'].boundary
-            map_boundary.plot(ax=ax, color=map_settings['color'])
-        
-        if 'farm' in kwargs:
-            self.gdf.loc[self.gdf['type']=='turbine'].plot(ax=ax, color=farm_settings['turbine']['color'], label='turbine')
-            self.gdf.loc[self.gdf['type']=='substation'].plot(ax=ax, color=farm_settings['oss']['color'], label='substation')
-        
-        ax.set_xlabel('Longitude')
-        ax.set_ylabel('Latitude')
-        ax.legend()
-
-        ax.set_xlim([-124.875, -124.55])
-        ax.set_ylim([40.025, 40.25])
-
-        fig.tight_layout()
-        
-        # Some GeoPandas Help
-        # to plot just one entry of a geoseries: gdf.loc[[0],'geometry'].plot()
-        # to get the columns of a gdf: gdf.columns
-        # merging gdf's
-        # adding columns to gdf's
-
-        return fig, ax
-
-
-
-    def addPoints(self, ax, pointlist=[], kwargs={}):
-
-        point_settings = kwargs['pointlist']
-
-        points = gpd.GeoDataFrame({'type':['nrel_channel','nrel_humboldt','nrel_crescent_city','hawaii'],
-                                   'geometry': pointlist}, crs='EPSG:4326')
-        
-        points.plot(ax=ax, color=point_settings['color'], marker=point_settings['marker'], label=point_settings['label'])
-    
-
-    def addState(self, ax, states=[], kwargs={}):
-        
-        for state in states:
-            state_settings=kwargs[state]
-
-            state_geom = self.usa.loc[self.usa['type']==state]
-            if 'boundary' in state_settings:
-                state_geom = state_geom.boundary
-
-            newstate = gpd.GeoDataFrame({'type':state, 'geometry':state_geom}, crs='EPSG:4326')
-
-            newstate.plot(ax=ax, color=state_settings['color'])
-
-
-
-
-    
 
-"""
 
 
 
