Solar Energy Mapping

A Walk Through of the Coding and Analysis to locate optimal Locations for a new Solar Energy Plant.

Introduction:

The research question this project aims to answer is where are the optimal locations for a new solar power plant in the state of Montana?. There are many factors that contribute to this, including solar irradiance potential, existing power infrastructure, areas that are in greater need of solar power, topography, and land use.

In stage one of this analysis, areas with and without access to power from one of the 8 existing power plants in Montana were aggregated by county using a spatial select with a service range polygon generated from analysis on a network dataset.

In this step in the analysis, this python notebook will be used to perform network optimizations to locate optimal areas for a new power plant that would service an area that does not yet have access to solar power.

To get started, lets take a brief look at the data sources used in this analysis.

Data Sources:

County Lines Layer

Montana Department of Commerce. (n.d.). Montana counties dataset. Montana State Library Open Data Portal. Retrieved October 16, 2024, from https://ceic-mtdoc.opendata.arcgis.com/datasets/d792b475378748b49bdb8023e555d45a_0/explore

Solar Irradiance Raster Layer

The World Bank Group. (n.d.). Global Solar Atlas – Solar data for the USA. Retrieved October 16, 2024, from https://globalsolaratlas.info/download/usa

Powerlines Polyline Layer

U.S. Energy Information Administration (EIA). (n.d.). U.S. electric transmission lines dataset. EIA Open Data Portal. Retrieved October 16, 2024, from https://atlas.eia.gov/datasets/bd24d1a282c54428b024988d32578e59_0/explore?location=45.751378%2C-111.100067%2C10.21

Existing Power Plants Layer

U.S. Energy Information Administration (EIA). (n.d.). Existing solar power plants in the USA dataset. EIA Open Data Portal. Retrieved October 16, 2024, from https://atlas.eia.gov/datasets/bf5c5110b1b944d299bb683cdbd02d2a_0/explore?location=45.454269%2C-108.987261%2C7.49

Energy Consumption by County

National Renewable Energy Laboratory (NREL). (2020). Net electricity and natural gas consumption by county (2020). SLOPE Data Viewer. Retrieved October 16, 2024, from https://maps.nrel.gov/slope/data-viewer?filters=%5B%5D&layer=energy-consumption.net-electricity-and-natural-gas-consumption&year=2020&res=county

Part One Graphic: Existing Infrastructure

Data Preparation:
import arcpyimport matplotlib.pyplot as pltimport numpy as np#read in gdb gdb_path = cwd + "Solar.gdb"
#set up vector layers #--------------------------------------------------------------------------------------------------------------------existing_solar_path = gdb_path + "\Power_Plants_Projected_Clipped" #power plant point layerexisting_solar = arcpy.MakeFeatureLayer_management(existing_solar_path, "existing_solar_layer")
powerlines_path = gdb_path + "\MontanaPowerlines_Clip" # powerline vector line layerpowerlines = arcpy.MakeFeatureLayer_management(powerlines_path, "Powerlines_layer")
census_path = gdb_path + "\CENSUS_2020_PL94171_MONTANA_COUNTY" #census county linescensus = arcpy.MakeFeatureLayer_management(census_path, "census_data_county")#This layer is in the wrong coordinate system (web mercator), it'll need to get reprojected
#reprojecting census data to correct Coordinate Systemcensus_reprojected_path = gdb_path + r"\Census_Reprojected"coord = arcpy.Describe(existing_solar).spatialReferencearcpy.Project_management( in_dataset= census_path, out_dataset = census_reprojected_path, out_coor_system = coord)#reassign projected data to the census variablecensus = arcpy.MakeFeatureLayer_management(census_reprojected_path, "census_reprojected")
counties_with_access_path = gdb_path + "\Counties_with_access_to_solar_power_from_an_existing_map" #counties with accessaccess = arcpy.MakeFeatureLayer_management(counties_with_access_path, "Counties_with_access")
service_area_path = gdb_path + "\Plant_Service_Area" #service area polygon from previous analysisservice = arcpy.MakeFeatureLayer_management(service_area_path, "Plant_service_area")

#set up raster layers #---------------------------------------------------------------------------------------------------------------------powerVoltageout_path = gdb_path + "\PVOUT_Projected_clipped"pvo = arcpy.Raster(powerVoltageout_path)
globalHorizontalIrradiation_path = gdb_path + "/GHI_Projected_Clipped"ghi = arcpy.Raster(globalHorizontalIrradiation_path)ghi.save(globalHorizontalIrradiation_path)

# Function to check spatial referencedef check_projection(layer): desc = arcpy.Describe(layer) return desc.spatialReference.name
# Check projectionsprint("Projection Data for all layers: \n")print(f"Existing Solar Plants: {check_projection(existing_solar)}")print(f"Powerlines: {check_projection(powerlines)}")print(f"Census Counties: {check_projection(census)}")print(f"Counties with Access: {check_projection(access)}")print(f"Service Areas: {check_projection(service)}")print(f"Raster PVOUT: {pvo.spatialReference.name}")print(f"Raster GHI: {ghi.spatialReference.name}")
#check raster informationprint("\n\nRaster Data Information: \n")print(f"{pvo.name} dimensions {pvo.width} x {pvo.height}")print(f"{ghi.name} dimensions {ghi.width} x {ghi.height}")

Note: The GHI Raster has a much higher resolution than the PVout Raster

Firstly, what is GHI and PVout? Global Horizontal Irradiance (GHI) is one of the most important metrics in assessing the estimated solar output of a site. It describes the total amount of solar irradiance that lands on a horizontal surface. This is the sum of two other metrics, Direct Normal Irradiation (DNI) and Diffuse Horizontal Irradiation (DIF), which measure direct sunlight and indirect sunlight scattered by the atmosphere, respectively. Given how DIF measures scattered light, a higher ratio of DIF/GHI indicates an area with more confounding materials, like cloud cover, water vapor content, and air pollution.

PVout is a metric given by the Global Solar Atlas that estimates Photovoltaiac output from a standard solar panel given the GHI, DNI, and DIF values for a region in kWh/kWp units. Here is a link to Global Solar Atlas' website that describes these metrics in detail : https://globalsolaratlas.info/support/faq

In this analysis, we will be using the GHI raster layer due to its significantly higher resolution. Also given how GHI is the key metric for PVout, this should yeild trustworthy results. We can use the PVout layer to help validate the results of this analysis even though this layer does not have a sufficent resolution for performing optimization.

Converting Suitability Raster to Polygon

Polygon data gathered from a suitability raster will be the most useful in determining optimal areas to place a new solar power plant because it will be easiest to combine with the other polygon and polyline features in this dataset.

Before we can convert the raster to a polygon, first we must convert the GHI layer from a float raster to an int raster. The arcpy library I am using for this analysis uses a function called RasterToPolygon_conversion(), which requires an input raster in int format.

This means the raster must be reclassified in integer format, without accidentally smudging the GHI values stored in float format.

Lets take a close look at the way the GHI float raster is currently classified.

# a function to Describe raster to get basic statsdef describe_stats(raster):    raster_desc = arcpy.Describe(raster)    raster_props = arcpy.sa.CellStatistics(raster.name, "MAXIMUM", "DATA")    #list basic stats    print(f"Raster Name: {raster_desc.name}")    print(f"Spatial Reference: {raster_desc.spatialReference.name}")    print(f"Min: {raster_props.minimum}")    print(f"Max: {raster_props.maximum}")    print(f"Mean: {raster_props.mean}")    print(f"Standard Deviation: {raster_props.standardDeviation}")

    # Convert raster to NumPy array to plot the stats    arr = arcpy.RasterToNumPyArray(raster, nodata_to_value=0) # set nodata values to 0, so integer raster case works        # Plot a histogram to visualize the classification scheme    plt.hist(arr[arr != 0].flatten(), bins=50, color='blue', alpha=0.7)    plt.title(f"{raster.name} Raster Value Distribution")    plt.xlabel(f"{raster.name} Value")    plt.ylabel("Frequency")
    #display the histogram.     plt.show()
# Generate report describing the stats from the GHI rasterdescribe_stats(ghi)

OUTPUTS:

Raster Name: GHI_Projected_Clipped

Spatial Reference: NAD_1983_UTM_Zone_12N

Min: 1.4509999752044678

Max: 4.581999778747559

Mean: 4.022209022653781

Standard Deviation: 0.2049334070089039

Now that we can see the classification scheme that is currently in use, a new classification in integer format can be dirived.

I've decided to reclassify the raster using quantiles as the data appears to be relatively normal. This will allow me to sort the raster values into five classes: low, medium low, medium, medium high, and high. Later on, only areas in medium high or high will be considered for an optimal site location.

I had many issues with the arcpy library functions for reclassification, as the Reclass() function requires an int raster in the first place. So, I wrote my own reclassification function

# a function to reclassify a float raster to integer format in a quantiles classificationdef rasterFloatToIntQuantiles(raster, reclass_range): ''' Parameters: raster: (a float raster object) , reclass_range: (a list of lists) the quantile bounds and new value for each quantile range Returns: the raster object the raster object's absolute filepath ''' #start with a null variable to allow starting with original raster reclass_raster = None #apply reclassification using loop for each quantile using the Con() Meathod for lower, upper, intvalue in reclass_range: if reclass_raster is None: #starting with original raster reclass_raster = Con((raster >= lower) & (raster <= upper), intvalue) else: # append new quantile reclassification on top of reclassed raster for each quantile reclass_raster = Con((raster >= lower) & (raster <= upper), intvalue, reclass_raster) # now save the reclassified raster name = f"{raster.name}_int" # Ensure no unwanted extensions path = f"{gdb_path}\\{name}" # Ensure it's saved in the geodatabase reclass_raster.save(path) return reclass_raster, path # return the new raster object and its path so it can be passed around outside this function
# get quantile breaks in a numpy array#using nan (not a number) pertcentiles to allow data in very high and very low percentiles to be includedquantiles = np.nanpercentile(arcpy.RasterToNumPyArray(ghi, nodata_to_value=0), [20, 40, 60, 80])
#create a reclassification scheme using these quantilesreclass_range = [ [float('-inf'), quantiles[0], 1], # low, assigned int value 1 (below 20th percentile) [quantiles[0], quantiles[1], 2], # Medium Low, assigned int value 2 (between 20th and 40th percentile) [quantiles[1], quantiles[2], 3], # Medium, assigned int value 3 (between 40th and 60th percentile) [quantiles[2], quantiles[3], 4], # Medium High, assigned int value 4 (between 60th and 80th percentile) [quantiles[3], float('inf'), 5] # High, assigned int value 5 (above 80th percentile)]
#perform reclassification and conversionbuffered_suitability_int, buffered_suitability_int_path = rasterFloatToIntQuantiles(buffered_suitability, reclass_range)#describe the reclassified statistics to ensure data was not mutated in conversiondescribe_stats(buffered_suitability_int)

Here is the resulting classification scheme:

Quick Overview of the data so far

The figure above displays counties that already have access in white, counties that do not yet have access in grey, and the buffered suitability GHI values on a red to green color scheme. You can see that most of the green areas on the suitiblity raster cover areas that already have access to solar energy, which makes a lot of sense! The most optimal areas in the state have already been leveraged. The North West corner of the state has very poor annual GHI ratings, which offer a solid explanation for their status of "no access". That being said, there are certainly a few areas with high GHI that currently do not have access to solar yet. Most notabily, Gallatin County, Powder River County, and Carter County. (southwest sliver and south east counties)

Now we can begin searching for optimal sites based on population density data.

Next Step: Optimizing sites by proximity to population dense areas

import pandas as pd
#reset at top of block to allow rerunning
if "POPTOT_c" in census.columns:    census = census.drop(columns= ["POPTOT_c"])
#filter suitability layerhigh_suitability = ghi_output_clipped[ghi_output_clipped['gridcode'].isin([4, 5])]
#add in population data census tablecensusTable = pd.DataFrame(    arcpy.da.TableToNumPyArray(        gdb_path + "\CensusTable", # table path        ["NAME", "POPTOT"], #field names for population totals    ))print(census.columns)census = census.merge(censusTable, on="NAME", suffixes= ("_c", "_table")) #join table to layer
#add in power consumption power_use_path = gdb_path + "\energy_consumption_expenditure_business_as_usual_county"power_use = pd.DataFrame(    arcpy.da.TableToNumPyArray(        power_use_path, #TablePath        ["County_Name", "Sector", "Consumption_MMBtu"], #Field Names        where_clause= '"State_Name" = \'Montana\'' # Definition query for montana    ))# peek at power use data
power_use.head(5)

As you can see, the data table containing power usage data has many lines per county. Next, I'll group the data by county, and by sector. Then I'll generate a sum for each sector in each county for energy consumption.

# Group by County and Sector and sum consumptiongrouped = ( power_use.groupby(["County_Name", "Sector"])["Consumption_MMBtu"] .sum() .reset_index())
# Pivot the table to have sectors as columnsenergy_summary = grouped.pivot(index="County_Name", columns="Sector", values="Consumption_MMBtu").fillna(0)
# Add a total column energy_summary["Total"] = energy_summary.sum(axis=1)
# Round all values to an integerenergy_summary = energy_summary.round().astype(int)
# Preview the resultenergy_summary.head()

Alright, Now we have a data table by county that shows the total energy consumption for all three sections. Lets use this data to determine what counties are in greatest need of supplemental energy.

Weighted Energy Demand Index

We will achieve this by creating a Weighted Energy Demand Index Field for each county and then using a spatial select between high GHI areas. I am choosing to give equal weight to 4 variables, consumption by commercial, industrial, and residential sectors, and total population. As you can tell, given how residential and population are both considered, this index is aimed at the benefit of people and residential energy.

# Reindex energy_summary to match censuscensus = census.set_index(census['NAME'])energy_summary = energy_summary.reindex(census.index)print(census.index)print(energy_summary.index)
# Weighted Energy Demand Indexcensus['Energy_Index'] = ( energy_summary['commercial'] * 0.25 + # 25% weight on commercial sector consumption energy_summary['industrial'] * 0.25 + # 25% weight on industrial sector consumption energy_summary['residential'] * 0.25+ # 25% weight on residential sector consumption census['POPTOT_table'] * 0.25 # 25% weight on the total population)
#There was a slight error here where somehow the census layer lost its correct spatialReference. Quick Reset.census = census.to_crs(26912)
# This column causes errors in the join, it came from the first join operation we did and we don't need it anymoreif 'index_right' in high_suitability.columns: high_suitability = high_suitability.drop(columns=['index_right']) # rename column to avoid errorcensus = census.rename(columns={'NAME': 'census_NAME'})
# spatial joinhigh_suitability = gpd.sjoin(high_suitability, census, how="left", predicate="intersects")

Using the weighted demand index, I was able to identify a few spots that are could be valid locations for a new solar power plant in the state of Montana.

The blue X's mark locations that serve a new population whereas the Purple X's mark the most optimal locations regardless of existing infrastructure.

Conclusion:

I have high confidence in the performance of my optimization model due to the similarity of my detected optimal locations to the locations of existing solar power plants (Comparison below). The two blue x's identify the new optimal sites to bring solar power to a broader population. These two locations are nearby Bozeman and Missoula, two of the three most populated cities in the state. These area's make sense due to the high weight on bulk energy consumption compared to PVout metrics. One of the key caveats with this graphic is that the x's mark the centriods of the county, not the particular pixel identified as the optimal location. I could have gone back to locate the exact pixel identified as optimal, however I believe that the visual on this map allows for an easy indicator to pick a location considering factors that were ignored by the weighted energy index, like land ownership and construction feasibility. All of these data layers can be exported to .shp shapefile format to use in interactive map areas, this would be very useful to a planning official looking to investigate possible sites.

I hope you enjoyed reading some of my analysis!

Home

Next: Remote Sensing

Page updated

Google Sites

Report abuse