Python 有一些非常酷的地理空间库,但最新的一个是 PyGMT,它是著名的通用映射工具的强大包装器。在继续学习 PyGMT 的过程中,我发现自己需要绘制一张地质图。经过多次试验和错误,我能够在 PyGMT 社区的帮助下让它工作。让我们来看看在 Windows 上的 PyGMT v0.4.1 中使用预先存在的 shapefile 重新创建地质图的一种方法:

此演示将使用USGS 网站上提供的以下文件:

1.俄勒冈地质图 shapefiles 的压缩文件夹

2.地质单元的颜色信息文件(点击链接后,右键单击“另存为”)

在关注此演示时,请注意所有代码块都旨在作为完整代码的一部分运行。本页底部提供了完整的代码。

设置文件夹

要按照编写的这个演示进行操作,您需要编辑脚本以更改main_dir变量以反映您的 User 文件夹名称。完成此操作后,您将需要在该目录中创建一些文件夹,特别是 Oregon_Geologic_Map_Demo 文件夹,并在其中创建 Methods、Data 和 Results 文件夹;在 Data 文件夹中创建一个 Conditioned_SHP 文件夹。您可以手动创建它们或从任何地方运行演示脚本,如果它们不存在,它们将为您创建它们。创建文件夹后,将解压缩的 shapefile 文件夹 (ORgeol_dd) 和颜色信息文件 (lithrgb.txt) 复制到 Data 文件夹。

所需目录树的示例如下:


<img altu003d"image" srcsetu003d"/_next/image?urlu003dhttps%3A%2F%2Fcdn.hackernoon.com%2Fimages%2F-sj6d26r2.png&wu003d640&qu003d75 1x, /_next/image?url u003dhttps%3A%2F%2Fcdn.hackernoon.com%2Fimages%2F-sj6d26r2.png&wu003d1080&qu003d75 2x" srcu003d"/_next/image?urlu003dhttps%3A%2F%2Fcdn.hackernoon.com%2Fimages %2F-sj6d26r2.png&wu003d1080&qu003d75" 解码u003d"async" data-nimgu003d"intrinsic" styleu003d"position:absolute;top:0;left:0;bottom:0;right:0;box-sizing:边框;填充:0;边框:无;边距:自动;显示:块;宽度:0;高度:0;最小宽度:100%;最大宽度:100%;最小高度:100%;最大-height:100%;object-fit:contain" classu003d"image" loadingu003d"lazy">

    # Creates the project folders
    def Create_Folders(self):
        import os

        # Main directory path
        self.main_dir = r'C:\Users\USER\Desktop\Oregon_Geologic_Map_Demo'

        # Path to the directory holding the project data files
        data_folder = os.path.join(self.main_dir, 'Data')

        # Path to the directory holding the project Python scripts
        methods_folder = os.path.join(self.main_dir, 'Methods')

        # Path to the directory holding the map generated by the scripts
        results_folder = os.path.join(self.main_dir, 'Results')

        # Path to the directory holding the conditioned shapefile data
        conditioned_shapefiles_folder = os.path.join(self.main_dir, 'Data', 'Conditioned_SHP')

        directories = [self.main_dir, data_folder, methods_folder, results_folder, conditioned_shapefiles_folder]

        # Iterates through the list of directories and creates them if they don't already exist
        for directory in directories:
            os.makedirs(directory, exist_ok = True)
        
        exit()

演示脚本创建的数据文件将保存到_Data_文件夹,俄勒冈州的地质图将保存到_Results_文件夹。 Methods 文件夹应包含演示脚本。

**调节 Shapefile **

需要对 shapefile(.shp 扩展名)进行调节,以准备制作调色板表(.cpt 扩展名)。调色板表需要地质单位名称采用它可以理解的格式,这不需要空格或特殊字符。 shapefile 和调色板表需要具有相同的地质单位名称,否则为地质单位着色将不起作用,并且由于调色板表是限制因素,我们需要调整 shapefile 以便单位名称将是同意。为简单起见,我遵循 USGS 提供的说明,使用其 ArcMap 样式表为地质图着色,该样式表指示按 ROCKTYPE1 列对单位进行着色。下一步是隔离 shapefile 中的 ROCKTYPE1 列值并调整它们:

    # List that holds the inital unit names
    unit_names_list = []
   
    # Conditions the shapefile
    def Condition_Shapefile(self):
        import os
        import geopandas as gpd  
        import pandas as pd
        import numpy as np
        
        # Shapefile of geologic unit polygons
        geo_polygons = os.path.join(self.main_dir, 'Data', 'ORgeol_dd', 'orgeol_poly_dd.shp')

        # Reads the shapefile into a DataFrame
        df_geo_polygons = gpd.read_file(geo_polygons, driver = 'SHP')

        # Creates a numpy array of the unique values in the ROCKTYPE1 column
        unit_names = df_geo_polygons['ROCKTYPE1'].unique()

        # List of unit names as they initally appear in the shapefile
        self.unit_names_list = list(unit_names)

        # Copy of unit names that is to be conditioned and subsituted for the original names
        conditioned_unit_names = list(unit_names)

        # Index of each character as they are read
        index = -1

        # Keys are what need to be replaced in words, and values are what they will be replaced with
        replacements = {' ':'_', '/':'_', '-':'_', '(':'_', ')':'_', 'é':'e', 'mã©lange':'melange'}

        # Iterates through the list of unique geologic unit names from the ROCKTYPE1 column and conditions them to the desired format
        for name in conditioned_unit_names:
            index += 1
            for old_value, new_value in replacements.items():
                if old_value in name:
                    conditioned_unit_names[index] = name.lower().replace(old_value, new_value)

        # Replaces the geologic unit names of the ROCKTYPE1 column in the dataframe with the conditioned names
        for name, conditioned_name in zip(unit_names, conditioned_unit_names):
            df_geo_polygons['ROCKTYPE1'] = df_geo_polygons['ROCKTYPE1'].replace(name, conditioned_name)

        # Save name for the conditioned shapefile
        geo_polygons_conditioned = os.path.join(self.main_dir,'Data', 'Conditioned_SHP', 'OR_geology_polygons_CONDITIONED.shp')

        # Saves the DataFrame as an ESRI shapefile
        df_geo_polygons.to_file(geo_polygons_conditioned, driver = 'ESRI Shapefile')

创建调色板表

既然带有多边形数据的 shapefile 已经设置好了,就需要制作一个调色板表来为多边形着色。以制表符描述的文本文件形式提供的颜色信息文件几乎可以用作调色板表,但需要进行一些调整。 “代码”列需要删除,标题也是如此,地质单元名称列需要位于 RGB 数据之前。此外,地质单位名称需要以与 shapefile 相同的方式进行调节,并且 RGB 值需要用正斜杠分隔:

    # Dictionary that holds the unit names from the cpt-like file prior to conditioning and the conditioned rgb colors
    cpt_data_dictionary = {}

    # Creates a color palette table file
    def Create_Color_Palette_Table(self):
        import os
        import pandas as pd
        import re

        # Path to table of geologic unit colors
        geo_colors = os.path.join(self.main_dir, 'Data', 'lithrgb.txt')

        # Reads the table to a DataFrame, ignoring the "code" column and skipping the column names
        df_geo_colors = pd.read_csv(geo_colors, sep ='\t', usecols = [1,2,3,4], skiprows = 1, header = None)

        # Moves the geologic unit names column to be first
        df_geo_colors.set_index(df_geo_colors.columns[-1], inplace = True)

        # Resets the index so that the new column order is respected
        df_geo_colors.reset_index(inplace = True)

        # Save name for the new cpt
        geo_cpt = os.path.join(self.main_dir, 'Data', 'geology_color_palette.cpt')

        # Writes the DataFrame as a CSV file so that the data can be conditioned
        df_geo_colors.to_csv(geo_cpt, sep = '\t', header = None, index = False)

        with open(geo_cpt, 'r', encoding = 'latin-1') as f:
            data = f.read()

        # Keys are what need to be replaced in words, and values are what they will be replaced with
        replacements = {' ':'', '/':'_', '-':'_', '(':'_', ')':'_', 'é':'e', 'mã©lange':'melange'}

        # Iterates over the replacments dictionary and replaces the desired characters in the text
        for old_value, new_value in replacements.items():
            data = data.replace(old_value, new_value)

        # Converts the text to lowercase
        data = data.lower()

        # Regular expression that finds tab-spaces between numbers
        pattern = re.compile(r'(?<=\d)(\t)(?=\d)')

        # Uses the regular expression to replace tab-spaces with "/"
        data = pattern.sub('/', data)

        with open(geo_cpt, 'w') as f:
            f.write(data)


        # DataFrame used to hold the unit names from the cpt-like file prior to conditioning and the conditioned rgb colors
        self.df_cpt_data = pd.DataFrame()

        # Adds only the first column (unit name column) to a column in the new dataframe titled 'geo_unit'
        self.df_cpt_data['geo_unit'] = df_geo_colors.iloc[:, 0]

        # Reads in the conditioned cpt file using tab delimiters and no column titles
        df_cpt_data = pd.read_csv(geo_cpt, sep ='\t', header = None, encoding = 'latin-1')

        # Adds only the second column (unit color column) to a column in the new dataframe titled 'color'
        self.df_cpt_data['color'] = df_cpt_data.iloc[:, 1]

        # Converts the values of the 'geo_unit' and 'color' columns into a dictionary with 'geo_unit' as the key and 'color' as the value
        self.cpt_data_dictionary = pd.Series(self.df_cpt_data.color.values, index=self.df_cpt_data.geo_unit).to_dict()

创建 Postscript 图例文件

为了给地质图的颜色提供背景,需要创建一个图例;这将通过创建一个 postscript 文件来完成:

    # Creates a legend postscript file
    def Create_Legend(self):
        import os
        import pandas as pd

        # File path to the geologic unit legend postscript file
        geo_legend = os.path.join(self.main_dir, 'Data', 'geologic_unit_legend.txt')

        # Index counter for determining position in unit_names_list
        index = -1

        # Iterates through the list of unit names and capitalizes the first letter of each entry
        for name in self.unit_names_list:
            index += 1
            self.unit_names_list[index] = name.capitalize()

        # Lists of the unit names and colors to be included in the legend
        selected_unit_names = []
        selected_colors = []

        # Iterates through the dictionary of unit names and colors, capitalizes the first letter of each unit name, and then matches it with the list of inital unit names.
        # Then those unit names and colors are added to respective lists
        for name, color in self.cpt_data_dictionary.items():
            cap_name = name.capitalize()
            if cap_name in self.unit_names_list:
                selected_unit_names.append(cap_name)
                selected_colors.append(color)

        # Creates a postscript file and writes some explainer text and the column format
        with open(geo_legend, 'w') as f:
            f.write(
                '# G is vertical gap, V is vertical line, N sets # of columns,\n'
                '# D draws horizontal line, H is header, L is column header,\n'
                '# S is symbol,\n' 
                '# format of: symbol position, symbol type,\n' 
                '# format of: symbol size, symbol color, symbol border thickness, text position, text\n'

                'N 2\n'
            )

        # Iterates through the unit names and colors, and adds the symbol+text lines to the existing postscript file
        with open(geo_legend, 'a') as f:
            for color, name in zip(selected_colors, selected_unit_names):
                f.write(
                    'S 0.1i r 0.1i {} 0.1p 0.20i {}\n'
                    '\n'
                    .format(color, name)
                )

绘制地质图

现在地质图终于可以绘制了:

    # Plots the map    
    def Plot_Map(self):
        import os
        import geopandas as gpd
        import pandas as pd
        import pygmt


        # Map save name
        save_name = os.path.join(self.main_dir, 'Results', 'Oregon_Geologic_Map_Demo.png')

        # Geologic unit polygons
        geo_unit_data = os.path.join(self.main_dir, 'Data', 'Conditioned_SHP', 'OR_geology_polygons_CONDITIONED.shp')

        # Geologic unit color palette
        geo_unit_color = os.path.join(self.main_dir, 'Data', 'geology_color_palette.cpt')

        geo_unit_legend = os.path.join(self.main_dir, 'Data', 'geologic_unit_legend.txt')

        # Extent defining the area of interest (Oregon) <min lon><max lon><min<lat><max lat>
        region = [-126.418246, -116.462397, 41.984295, 46.297196]

        # Map projection (Mercator): <type><size><units>
        projection = 'M6i' 

        # Frame  annotations: [<frame sides with/without ticks>, <x-axis tick rate>, <y-axis tick rate>]
        frame = ['SWne', 'xa', 'ya'] 

        # Polygon outline pens: <size><color>
        pens = {'geology':'0.1p,black'} 

        # Transparency of layer, transparency of save file background
        transparency = {'geology':50, 'save_file':False}    

        df_geo_polygons = gpd.read_file(geo_unit_data, driver = 'SHP')

        # Establishes figure to hold map layers
        fig = pygmt.Figure()

        # Forces ticks to display as degree decimal
        pygmt.config(FORMAT_GEO_MAP = 'D')

        # Forces the map frame annotation to be smaller
        pygmt.config(FONT_ANNOT_PRIMARY = '8p,Helvetica,black')


        # Basemap layer
        fig.basemap(
            region = region,
            projection = projection,
            frame = frame
            )

        # Geologic unit layer
        fig.plot(
            # File data - automatically detects polygon coordinates if in last column
            data = df_geo_polygons,
            # Sets polygon outline colour  
            pen = pens['geology'],
            # Sets polygon color map
            cmap = geo_unit_color,
            # Sets color to vary with selected column
            color = '+z',
            # Force close polygons
            close = True,
            # Sets the column used to map polygon colors (in this case colors polygons by name of geologic unit). Column name appears to be lowercase as a product of conditioning
            aspatial = 'Z=ROCKTYPE1',
            # Sets layer transperancy
            transparency = transparency['geology'],
            # Commandline feedback for debugging 
            #verbose=True,  
            )

        # Plots the coastlines and political boundaries
        fig.coast(
            # Displays national boundaries (1) with 0.8 point gray40 lines, and does the same for state boundaries (2)
            borders = ['1/0.8p,gray40', '2/0.8p,gray40'], 
            # Displays coast outlines in 0.3 point black lines, and lakeshore outlines in .1 point black lines
            shorelines = ['1/0.3p,black', "2/0.1p,black"],
            # Sets resolution full (f) [highest setting]    
            resolution = 'f',
            # Sets water color
            water = 'lightskyblue2',  
            )

        # Plots a legend of the geologic unit names and respective colors
        fig.legend(
            spec = geo_unit_legend, # pslegend file
            position = 'jBL+o15.5/-4c+w10/12c', # plots text justifed bottom left (jBL) and offsets (+o) it by 15.5cm on the x-axis and -4cm on the y-axis (15.5/-4c), and establises width of columns(?)/legend area(?) (+w) as 10cm on the x-axis and 12cm on the y-axis (10/12c)
        )

        # Saves a copy of the generated figure
        fig.savefig(save_name, transparent = transparency['save_file'])

结果

瞧,俄勒冈州的地质图是通过 PyGMT 用 Python 渲染的:


<img altu003d"源数据:Walker, G.W. 和 MacLeod, N.S.,1991 年,俄勒冈州地质图:美国地质调查局,比例 1:500,000" srcsetu003d"/_next/image?urlu003dhttps%3A%2F%2Fcdn .hackernoon.com%2Fimages%2F-cibb26rt.png&wu003d828&qu003d75 1x, /_next/image?urlu003dhttps%3A%2F%2Fcdn.hackernoon.com%2Fimages%2F-cibb26rt.png&wu003d1920&qu003d75 2x " srcu003d"/_next/image?urlu003dhttps%3A%2F%2Fcdn.hackernoon.com%2Fimages%2F-cibb26rt.png&wu003d1920&qu003d75" 解码u003d"async" data-nimgu003d"intrinsic" styleu003d "位置:绝对;顶部:0;左侧:0;底部:0;右侧:0;框尺寸:边框框;填充:0;边框:无;边距:自动;显示:块;宽度:0;高度:0;最小宽度:100%;最大宽度:100%;最小高度:100%;最大高度:100%;对象适配:包含" classu003d"image" loadingu003d"lazy">

源数据:Walker, G.W.和 MacLeod, N.S.,1991 年,俄勒冈州地质图:美国地质调查局,比例 1:500,000

可以采取的另一个步骤是使用 DEM 支撑地图以提供更多背景信息。但是,这只是稍微超出了本演示的预期范围。我很高兴学习如何将公开可用的地质图转换为我可以在 Python 中使用的形式,我希望这可以帮助其他希望做同样事情的人。

完整代码

'''
This script builds a geologic map of Oregon using an existing shapefile

Data source website: https://mrdata.usgs.gov/geology/state/state.php?state=OR
Shapefile source: http://pubs.usgs.gov/of/2005/1305/data/ORgeol_dd.zip
Color info file (right-click and "save page as"): https://mrdata.usgs.gov/catalog/lithrgb.txt

Setup instructions:
1) To run this script as-is, you will need to change the main_dir variable to reflect your desired path. 
2) You will need to create some folders at the destination, specifically a "Oregon_Geologic_Map_Demo" folder, and within it the "Methods", "Data", and "Results" folders.
You can either do that manually or run this script once from anywhere to automatically create them. 
3) Once you have created them, move the unzipped shapefile folder titled ORgeol_dd to the Data folder, as well as the color info file.
4) Now, save this script to the Methods folder and run it from there, if you wish.

Data files that are generated by this script will save to the Data folder, and the resulting map will save to Results
'''


# Controls whether the project folders data are automatically created
create_folders = False

# Controls whether the shapefile data is conditioned
condition_shp = True
    
# Controls whether a color palette table is made
create_cpt = True

# Controls whether a geologic unit legend postscript file is created
create_legend = True

# Controls whether the geologic map is created
plot_map = True


# Class that holds all the functions pretaining to creating a geologic map of Oregon
class Map_Maker():
    
    # Main directory path
    main_dir = r'C:\Users\USER\Desktop\Oregon_Geologic_Map_Demo'

    # Creates the project folders
    def Create_Folders(self):
        import os

        # Path to the directory holding the project data files
        data_folder = os.path.join(self.main_dir, 'Data')

        # Path to the directory holding the project Python scripts
        methods_folder = os.path.join(self.main_dir, 'Methods')

        # Path to the directory holding the map generated by the scripts
        results_folder = os.path.join(self.main_dir, 'Results')

        # Path to the directory holding the conditioned shapefile data
        conditioned_shapefiles_folder = os.path.join(self.main_dir, 'Data', 'Conditioned_SHP')

        directories = [self.main_dir, data_folder, methods_folder, results_folder, conditioned_shapefiles_folder]

        # Iterates through the list of directories and creates them if they don't already exist
        for directory in directories:
            os.makedirs(directory, exist_ok = True)
        
        exit()

    

    # List that holds the inital unit names
    unit_names_list = []

    # Conditions the shapefile
    def Condition_Shapefile(self):
        import os
        import geopandas as gpd  
        import pandas as pd
        import numpy as np
        
        # Shapefile of geologic unit polygons
        geo_polygons = os.path.join(self.main_dir, 'Data', 'ORgeol_dd', 'orgeol_poly_dd.shp')

        # Reads the shapefile into a DataFrame
        df_geo_polygons = gpd.read_file(geo_polygons, driver = 'SHP')

        # Creates a numpy array of the unique values in the ROCKTYPE1 column
        unit_names = df_geo_polygons['ROCKTYPE1'].unique()

        # List of unit names as they initally appear in the shapefile
        self.unit_names_list = list(unit_names)

        # Copy of unit names that is to be conditioned and subsituted for the original names
        conditioned_unit_names = list(unit_names)

        # Index of each character as they are read
        index = -1

        # Keys are what need to be replaced in words, and values are what they will be replaced with
        replacements = {' ':'_', '/':'_', '-':'_', '(':'_', ')':'_', 'é':'e', 'mã©lange':'melange'}

        # Iterates through the list of unique geologic unit names from the ROCKTYPE1 column and conditions them to the desired format
        for name in conditioned_unit_names:
            index += 1
            for old_value, new_value in replacements.items():
                if old_value in name:
                    conditioned_unit_names[index] = name.lower().replace(old_value, new_value)

        # Replaces the geologic unit names of the ROCKTYPE1 column in the dataframe with the conditioned names
        for name, conditioned_name in zip(unit_names, conditioned_unit_names):
            df_geo_polygons['ROCKTYPE1'] = df_geo_polygons['ROCKTYPE1'].replace(name, conditioned_name)

        # Save name for the conditioned shapefile
        geo_polygons_conditioned = os.path.join(self.main_dir,'Data', 'Conditioned_SHP', 'OR_geology_polygons_CONDITIONED.shp')

        # Saves the DataFrame as an ESRI shapefile
        df_geo_polygons.to_file(geo_polygons_conditioned, driver = 'ESRI Shapefile')



    # Dictionary that holds the unit names from the cpt-like file prior to conditioning and the conditioned rgb colors
    cpt_data_dictionary = {}

    # Creates a color palette table file
    def Create_Color_Palette_Table(self):
        import os
        import pandas as pd
        import re

        # Path to table of geologic unit colors
        geo_colors = os.path.join(self.main_dir, 'Data', 'lithrgb.txt')

        # Reads the table to a DataFrame, ignoring the "code" column and skipping the column names
        df_geo_colors = pd.read_csv(geo_colors, sep ='\t', usecols = [1,2,3,4], skiprows = 1, header = None)

        # Moves the geologic unit names column to be first
        df_geo_colors.set_index(df_geo_colors.columns[-1], inplace = True)

        # Resets the index so that the new column order is respected
        df_geo_colors.reset_index(inplace = True)

        # Save name for the new cpt
        geo_cpt = os.path.join(self.main_dir, 'Data', 'geology_color_palette.cpt')

        # Writes the DataFrame as a CSV file so that the data can be conditioned
        df_geo_colors.to_csv(geo_cpt, sep = '\t', header = None, index = False)

        with open(geo_cpt, 'r', encoding = 'latin-1') as f:
            data = f.read()

        # Keys are what need to be replaced in words, and values are what they will be replaced with
        replacements = {' ':'', '/':'_', '-':'_', '(':'_', ')':'_', 'é':'e', 'mã©lange':'melange'}

        # Iterates over the replacments dictionary and replaces the desired characters in the text
        for old_value, new_value in replacements.items():
            data = data.replace(old_value, new_value)

        # Converts the text to lowercase
        data = data.lower()

        # Regular expression that finds tab-spaces between numbers
        pattern = re.compile(r'(?<=\d)(\t)(?=\d)')

        # Uses the regular expression to replace tab-spaces with "/"
        data = pattern.sub('/', data)

        with open(geo_cpt, 'w') as f:
            f.write(data)


        # DataFrame used to hold the unit names from the cpt-like file prior to conditioning and the conditioned rgb colors
        self.df_cpt_data = pd.DataFrame()

        # Adds only the first column (unit name column) to a column in the new dataframe titled 'geo_unit'
        self.df_cpt_data['geo_unit'] = df_geo_colors.iloc[:, 0]

        # Reads in the conditioned cpt file using tab delimiters and no column titles
        df_cpt_data = pd.read_csv(geo_cpt, sep ='\t', header = None, encoding = 'latin-1')

        # Adds only the second column (unit color column) to a column in the new dataframe titled 'color'
        self.df_cpt_data['color'] = df_cpt_data.iloc[:, 1]

        # Converts the values of the 'geo_unit' and 'color' columns into a dictionary with 'geo_unit' as the key and 'color' as the value
        self.cpt_data_dictionary = pd.Series(self.df_cpt_data.color.values, index=self.df_cpt_data.geo_unit).to_dict()



    # Creates a legend postscript file
    def Create_Legend(self):
        import os
        import pandas as pd

        # File path to the geologic unit legend postscript file
        geo_legend = os.path.join(self.main_dir, 'Data', 'geologic_unit_legend.txt')

        # Index counter for determining position in unit_names_list
        index = -1

        # Iterates through the list of unit names and capitalizes the first letter of each entry
        for name in self.unit_names_list:
            index += 1
            self.unit_names_list[index] = name.capitalize()

        # Lists of the unit names and colors to be included in the legend
        selected_unit_names = []
        selected_colors = []

        # Iterates through the dictionary of unit names and colors, capitalizes the first letter of each unit name, and then matches it with the list of inital unit names.
        # Then those unit names and colors are added to respective lists
        for name, color in self.cpt_data_dictionary.items():
            cap_name = name.capitalize()
            if cap_name in self.unit_names_list:
                selected_unit_names.append(cap_name)
                selected_colors.append(color)

        # Creates a postscript file and writes some explainer text and the column format
        with open(geo_legend, 'w') as f:
            f.write(
                '# G is vertical gap, V is vertical line, N sets # of columns,\n'
                '# D draws horizontal line, H is header, L is column header,\n'
                '# S is symbol,\n' 
                '# format of: symbol position, symbol type,\n' 
                '# format of: symbol size, symbol color, symbol border thickness, text position, text\n'

                'N 2\n'
            )

        # Iterates through the unit names and colors, and adds the symbol+text lines to the existing postscript file
        with open(geo_legend, 'a') as f:
            for color, name in zip(selected_colors, selected_unit_names):
                f.write(
                    'S 0.1i r 0.1i {} 0.1p 0.20i {}\n'
                    '\n'
                    .format(color, name)
                )



    # Plots the map    
    def Plot_Map(self):
        import os
        import geopandas as gpd
        import pandas as pd
        import pygmt


        # Map save name
        save_name = os.path.join(self.main_dir, 'Results', 'Oregon_Geologic_Map_Demo.png')

        # Geologic unit polygons
        geo_unit_data = os.path.join(self.main_dir, 'Data', 'Conditioned_SHP', 'OR_geology_polygons_CONDITIONED.shp')

        # Geologic unit color palette
        geo_unit_color = os.path.join(self.main_dir, 'Data', 'geology_color_palette.cpt')

        geo_unit_legend = os.path.join(self.main_dir, 'Data', 'geologic_unit_legend.txt')

        # Extent defining the area of interest (Oregon) <min lon><max lon><min<lat><max lat>
        region = [-126.418246, -116.462397, 41.984295, 46.297196]

        # Map projection (Mercator): <type><size><units>
        projection = 'M6i' 

        # Frame  annotations: [<frame sides with/without ticks>, <x-axis tick rate>, <y-axis tick rate>]
        frame = ['SWne', 'xa', 'ya'] 

        # Polygon outline pens: <size><color>
        pens = {'geology':'0.1p,black'} 

        # Transparency of layer, transparency of save file background
        transparency = {'geology':50, 'save_file':False}    

        df_geo_polygons = gpd.read_file(geo_unit_data, driver = 'SHP')

        # Establishes figure to hold map layers
        fig = pygmt.Figure()

        # Forces ticks to display as degree decimal
        pygmt.config(FORMAT_GEO_MAP = 'D')

        # Forces the map frame annotation to be smaller
        pygmt.config(FONT_ANNOT_PRIMARY = '8p,Helvetica,black')


        #  Basemap layer
        fig.basemap(
            region = region,
            projection = projection,
            frame = frame
            )

        # Geologic unit layer
        fig.plot(
            # data file - automatically detects polygon coordinates if in last column
            data = df_geo_polygons,
            # Sets polygon outline colour  
            pen = pens['geology'],
            # Sets polygon color map
            cmap = geo_unit_color,
            # Sets color to vary with selected column
            color = '+z',
            # Force close polygons
            close = True,
            # Sets the column used to map polygon colors (in this case colors polygons by name of geologic unit). Cloumn name appears to be lowercase as a product of conditioning
            aspatial = 'Z=ROCKTYPE1',
            # Sets layer transperancy
            transparency = transparency['geology'],
            # Commandline feedback for debugging 
            #verbose=True,  
            )

        # Plots the coastlines and political boundaries
        fig.coast(
            # Displays national boundaries (1) with 0.8 point gray40 lines, and does the same for state boundaries (2)
            borders = ['1/0.8p,gray40', '2/0.8p,gray40'], 
            # Displays coast outlines in 0.3 point black lines, and lakeshore outlines in .1 point black lines
            shorelines = ['1/0.3p,black', "2/0.1p,black"],
            # Sets resolution full (f) [highest setting]    
            resolution = 'f',
            # Sets water color
            water = 'lightskyblue2',  
            )

        # Plots a legend of the geologic unit names and respective colors
        fig.legend(
            spec = geo_unit_legend, # pslegend file
            position = 'jBL+o15.5/-4c+w10/12c', # plots text justifed bottom left (jBL) and offsets (+o) it by 15.5cm on the x-axis and -4cm on the y-axis (15.5/-4c), and establises width of columns(?)/legend area(?) (+w) as 10cm on the x-axis and 12cm on the y-axis (10/12c)
        )

        # Saves a copy of the generated figure
        fig.savefig(save_name, transparent = transparency['save_file'])



data = Map_Maker()

function_controls = [create_folders, condition_shp, create_cpt, create_legend, plot_map]
functions = [data.Create_Folders, data.Condition_Shapefile, data.Create_Color_Palette_Table, data.Create_Legend, data.Plot_Map]

for control, function in zip(function_controls, functions):
    if control == True:
        function()