datasets
Every data source used in this course, and why it was chosen
This page documents every dataset the Solar Site Suitability project relies on: where it comes from, what it actually represents, why it was chosen over alternatives, and how it is downloaded. The code shown here is the exact acquisition script used to prepare the course data before Session 01. It is shown for reference and is not executed on this page, since running it requires personal credentials (a Google Earth Engine project) and produces large files that do not belong on a documentation page.
You will need a free Google Earth Engine account and a Google Cloud project ID, along with the Python packages earthengine-api, geemap, geopandas, osmnx, shapely, and python-dotenv.
Study Area Boundary
Every other dataset in this project is downloaded relative to a single boundary polygon. This boundary is generated by geocoding a place name through OpenStreetMap’s Nominatim service, which OSMnx wraps in a single function call. Geocoding a place name is the fastest way to get started, but it returns the official boundary OpenStreetMap has on record for that place, so it is always worth plotting the result and checking its size before relying on it for the rest of the project.
from pathlib import Path
import geopandas as gpd
import osmnx as ox
DATA_DIR = Path("data")
DATA_DIR.mkdir(parents=True, exist_ok=True)
PLACE_NAME = "Desierto de Tabernas, Almeria, Spain"
boundary_gdf = ox.geocode_to_gdf(PLACE_NAME)
boundary_gdf.plot()boundary_gdf.to_file(DATA_DIR / "boundary.gpkg", driver="GPKG")
polygon = boundary_gdf.geometry.iloc[0]ox.geocode_to_gdf already returns the boundary in EPSG:4326 (plain latitude and longitude), which is exactly the coordinate system Earth Engine expects when a region is passed to it. Reprojecting into a metric CRS happens later, in Session 03 and Session 07, once we actually need to measure distances or areas in meters.
Digital Elevation Model, Copernicus GLO-30
The DEM provides elevation values across the study area and is the source for two of the four suitability factors used later in the project: slope and aspect. This course uses Copernicus GLO-30, a global elevation model at 30 meter resolution, produced from radar data collected during the TanDEM-X mission between 2010 and 2015. It is one of the most widely used free global DEMs available today, alongside alternatives such as SRTM and ALOS World 3D, and is accessible directly through Google Earth Engine under the collection id COPERNICUS/DEM/GLO30_2024_1.
The original Earth Engine collection COPERNICUS/DEM/GLO30 has been deprecated and replaced by COPERNICUS/DEM/GLO30_2024_1. The year in the new name refers to when this edition of the dataset was reprocessed and republished, not to when the elevation was measured. The underlying radar data is still from 2010 to 2015. For natural terrain that has not been reshaped by mining, landslides, or large construction, this age difference has no meaningful effect on slope or aspect calculations. More detail is available on the Copernicus DEM overview page.
Technically, GLO30 is a Digital Surface Model rather than a Digital Terrain Model, meaning that in areas with buildings or dense vegetation, the elevation includes the height of what sits on top of the ground, not only the bare earth. Since this project specifically targets barren or low vegetation land for solar panel placement, this distinction rarely causes a practical problem.
import geopandas as gpd
import ee
import geemap
ee.Authenticate()
ee.Initialize(project="Your-ProjectID")
aoi = gpd.read_file(DATA_DIR / "boundary.gpkg")
aoi = geemap.geopandas_to_ee(aoi).geometry()dem_collection = ee.ImageCollection("COPERNICUS/DEM/GLO30_2024_1")
dem = dem_collection.select("DEM").mosaic().clip(aoi)
dem_output_path = DATA_DIR / "dem.tif"
geemap.ee_export_image(
dem,
filename=str(dem_output_path),
scale=30,
region=aoi,
file_per_band=False,
)geemap.ee_export_image downloads the raster straight to local disk and works well for small study areas, but Earth Engine limits how large a single direct download can be, roughly a few hundred megabytes. For a study area covering an entire province rather than a few square kilometers, this limit is exceeded and the export must instead go through ee.batch.Export.image.toDrive, which has no practical size limit but requires retrieving the finished file from Google Drive afterward, and often splits large exports into several tiles that need to be merged back together with rasterio.merge. The Sentinel-2 section below shows both approaches side by side.
Sentinel-2 Surface Reflectance Imagery
Sentinel-2 provides the multispectral imagery used to calculate NDVI, which distinguishes vegetated land from barren or sparsely vegetated land, one of the four suitability factors in this project. This course uses the Sentinel-2 Level-2A collection, which is already corrected for atmospheric effects (surface reflectance rather than top-of-atmosphere reflectance), avoiding the need to apply that correction manually. On Earth Engine this is the COPERNICUS/S2_SR_HARMONIZED collection.
A single Sentinel-2 scene is rarely usable on its own because of cloud cover, so the acquisition script filters the full image collection down to scenes with less than a chosen cloud cover percentage over a date range, masks out any remaining cloud and cloud shadow pixels using the Scene Classification Layer (SCL) band, and then takes the pixel-wise median across all remaining scenes. The result is a single, mostly cloud-free composite image rather than one individual acquisition, which is far more reliable for a fixed teaching dataset than trying to find one perfect scene by hand.
def mask_s2_clouds(image):
scl = image.select("SCL")
mask = scl.neq(3).And(scl.neq(8)).And(scl.neq(9)).And(scl.neq(10))
return image.updateMask(mask)
DATE_START = "2024-01-01"
DATE_END = "2024-12-31"
MAX_CLOUD_PERCENT = 10
s2_collection = (
ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
.filterBounds(aoi)
.filterDate(DATE_START, DATE_END)
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", MAX_CLOUD_PERCENT))
.map(mask_s2_clouds)
)
print(f"Number of matching images: {s2_collection.size().getInfo()}")s2_image = s2_collection.median().clip(aoi)
s2_export = s2_image.select(["B2", "B3", "B4", "B8"])Only four of the twelve available Sentinel-2 bands are exported: B2 (blue), B3 (green), B4 (red), and B8 (near infrared). Red and near infrared are the two bands NDVI actually needs, while blue and green are kept along for true color visualization purposes in later sessions. A full band description is available on the official Sentinel-2 bands reference.
The Scene Classification Layer assigns every pixel a class code. This project excludes class 3 (cloud shadow), class 8 (cloud, medium probability), class 9 (cloud, high probability), and class 10 (thin cirrus). Pixels in any of these classes are set to no data rather than included in the composite, since they do not represent the true reflectance of the ground.
sentinel_output_path = DATA_DIR / "sentinel.tif"
geemap.ee_export_image(
s2_export,
filename=str(sentinel_output_path),
scale=10,
region=aoi,
crs="EPSG:4326",
file_per_band=False,
)For a study area large enough to exceed the direct download limit described earlier, the same export instead goes through Google Drive:
task = ee.batch.Export.image.toDrive(
image=s2_export,
description="Sentinel",
folder="data_out",
region=aoi,
scale=10,
crs="EPSG:4326",
maxPixels=1e13,
)
task.start()task.start() only queues the export; it does not block or return the file. Progress has to be checked from the Earth Engine Code Editor’s Tasks tab or by polling task.status(), and the finished file then has to be downloaded from Google Drive manually or with the Google Drive API.
Road Network, OpenStreetMap
Roads represent one half of the infrastructure accessibility factor: land close to existing roads is more practical and less costly to develop for solar installations, since it is easier to reach for construction and maintenance. This course retrieves roads through OSMnx, a Python library that queries OpenStreetMap data through the Overpass API and returns it directly as a GeoDataFrame, without needing to manually download and clip an extract of an entire country.
road_graph = ox.graph_from_polygon(polygon, network_type="drive")
_, edges = ox.graph_to_gdfs(road_graph)
edges = gpd.clip(edges, boundary_gdf)ax = boundary_gdf.plot(facecolor="None", edgecolor="red")
edges.plot(ax=ax)edges.to_file(DATA_DIR / "roads.gpkg", driver="GPKG")The network_type="drive" argument restricts the result to roads that vehicles can travel on, excluding footpaths, cycleways, and similar minor ways that are irrelevant to construction access. OpenStreetMap coverage quality varies significantly by region: countries with an active mapping community, much of Western Europe among them, tend to have dense and reliable road data, while more remote areas may have gaps. It is worth visually inspecting the downloaded roads before relying on them for analysis, exactly as the plot above does.
OSMnx queries a public Overpass API server by default, which occasionally becomes temporarily unreachable under heavy load. If a connection error appears, switching to a mirror server, for example by setting ox.settings.overpass_endpoint = "https://overpass.kumi.systems/api/interpreter" before making the request, usually resolves it. More endpoints are listed on the Overpass API wiki page.
Power Line Network, OpenStreetMap
Power infrastructure represents the second half of the infrastructure accessibility factor: proximity to the existing electrical grid directly affects how expensive it is to connect a new solar installation. This data also comes from OpenStreetMap, using OSMnx’s more general feature query function rather than the road-specific graph function, since power infrastructure is not a routable network but a collection of tagged features. Two separate queries are made, one for power lines and one for the poles and towers that carry them, since both are useful for a complete picture of the grid, and each has a different geometry type.
power_tags = {"power": ["line", "minor_line"]}
power_line = ox.features_from_polygon(polygon, tags=power_tags)
power_line = gpd.clip(power_line, boundary_gdf)fig, ax = plt.subplots(figsize=(12, 7))
boundary_gdf.plot(ax=ax, facecolor="None", edgecolor="red")
edges.plot(ax=ax, label="Road")
power_line.plot(ax=ax, color="yellow", label="Power Line")
ax.legend()
plt.tight_layout()OpenStreetMap tags high voltage transmission lines as power=line and lower voltage local distribution lines as power=minor_line. Both are retrieved here since either type of line represents a usable grid connection point for this project. Full tagging conventions for power infrastructure are documented on the OpenStreetMap power wiki page.
tower_tags = {"power": ["tower", "pole"]}
towers = ox.features_from_polygon(polygon, tags=tower_tags)fig, ax = plt.subplots(figsize=(12, 7))
boundary_gdf.plot(ax=ax, facecolor="None", edgecolor="red")
edges.plot(ax=ax, label="Road")
power_line.plot(ax=ax, color="yellow", label="Power Line")
towers.plot(ax=ax, color="green", label="Towers")
ax.legend()
plt.tight_layout()Towers and poles are point features, so they cannot share a geometry column with the line features above. Rather than discarding them, both feature types are saved into the same GeoPackage as two separate layers, lines and towers, keeping everything related to power infrastructure in a single file.
power_line.to_file(DATA_DIR / "power_lines.gpkg", driver="GPKG", layer="lines")
towers.to_file(DATA_DIR / "power_lines.gpkg", driver="GPKG", layer="towers")Since power_lines.gpkg now holds two layers, reading it with a plain gpd.read_file("data/power_lines.gpkg") only returns the first one. The layer argument has to be given explicitly, for example gpd.read_file("data/power_lines.gpkg", layer="towers"), exactly the same pattern used for power.gpkg starting in Session 02.
Summary of Files Produced
Running the full acquisition script produces the following files inside the project’s data folder, all sharing the same study area boundary and coordinate reference system, ready to be used starting in Session 01.
| File | Source | Format | Resolution |
|---|---|---|---|
boundary.gpkg |
OSM geocoding | GeoPackage | Not applicable |
dem.tif |
Copernicus GLO-30 via Earth Engine | GeoTIFF | 30 m |
sentinel.tif |
Sentinel-2 L2A via Earth Engine | GeoTIFF | 10 m |
roads.gpkg |
OpenStreetMap via OSMnx | GeoPackage | Vector |
power_lines.gpkg (layer lines) |
OpenStreetMap via OSMnx | GeoPackage | Vector |
power_lines.gpkg (layer towers) |
OpenStreetMap via OSMnx | GeoPackage | Vector |
The size of the boundary polygon directly determines file size and download time for every dataset above. A boundary covering an entire province can produce a Sentinel-2 export of a gigabyte or more, split across several tiles, while a boundary a few kilometers across produces files of only a few megabytes and downloads in seconds. Keeping the boundary small is not just a matter of convenience, it keeps the whole course repository lightweight and every processing step in later sessions fast enough to run live in class.
Going Beyond This Course’s Datasets
The datasets described on this page were deliberately kept simple: a single free elevation model, one cloud-filtered Sentinel-2 composite, and OpenStreetMap layers for roads and power infrastructure. This choice is intentional. The course is meant to be approachable regardless of your background, so every dataset here can be downloaded in minutes, without specialized access requests, licensing negotiations, or heavy preprocessing. None of it requires prior remote sensing expertise to obtain or interpret.
That said, this is far from the only reasonable way to build a solar site suitability model. The academic and industry literature on renewable energy site selection uses a considerably richer set of inputs, and if you want to push this project further, it is worth exploring some of them. A few directions to look into:
- Solar irradiance data, such as the Global Solar Atlas or Copernicus Atmosphere Monitoring Service products, which measure actual sunlight availability rather than assuming it is uniform across the study area.
- Higher resolution or radar-derived DEMs, such as national LiDAR-based elevation datasets where available, which can resolve terrain features far more precisely than a 30 meter global model.
- Land cover and land use datasets, such as ESA WorldCover or Corine Land Cover, which classify land more precisely than an NDVI threshold and can rule out protected areas, agricultural land, or urban zones directly.
- Grid capacity and substation data, which national grid operators sometimes publish, going beyond simple proximity to a power line to also account for whether that line actually has spare capacity.
- Multi-criteria decision analysis methods, such as the Analytic Hierarchy Process (AHP), which are commonly used in the literature to justify the weights assigned to each suitability factor, rather than choosing them by hand as this course does.
A good starting point is searching for terms such as “solar site suitability GIS AHP” or “renewable energy site selection multi-criteria” in Google Scholar, which surfaces a large number of peer-reviewed studies describing exactly which datasets and weighting schemes they used for similar projects in different countries.