This guide has 2 routes. Pick one, or read both. Each route is complete on its own.
Sentinel-2 is a pair of European satellites that image most of the world's land at 10 m resolution, and the data is free. Getting it is the first step of almost every remote sensing task, and it is where beginners lose the most time: the wrong product, a cloudy scene, a 1 GB zip when ten megabytes would have done.
This guide shows two routes to the same data. Route A uses the Copernicus Browser, a website where you click. Route B uses a short Python script. Use route A when you want one or two whole scenes and like clicking. Use route B when you want a small area, many dates, or a result you can repeat.
Decide what to download
Three choices decide what you get.
The processing level. Level-1C is top-of-atmosphere reflectance, which still contains the haze between the satellite and the ground. Level-2A is surface reflectance, corrected for the atmosphere, and it comes with a scene classification layer (SCL) that marks cloud, shadow, water and more. For vegetation, water and land cover work, choose Level-2A. Level-1C is for people who want to run their own atmospheric correction.
The tile. Sentinel-2 data is cut into square tiles of about 110 km. A tile name looks like T43REQ. The 43 is the UTM zone, the R is a band of latitude, and EQ is a 100 km square inside that zone. One study area can fall in one tile or in two.
The date and the cloud. Clouds hide the ground, so you want a scene with little cloud over your area. A scene can be marked 30 percent cloudy over the whole tile and still be clear over your small corner, which is why you always look at the preview.
Route A: the Copernicus Browser
The Copernicus Data Space Ecosystem runs a free browser for all Sentinel data. The official documentation is linked at the end of this guide.
- Open
browser.dataspace.copernicus.eu. A welcome window offers Log in or Use anonymously. The official documentation says you need a free account to use the Browser, and downloading needs one. Register for free with Log in, then sign in. - Find your place. The box at the top right of the map, Go to Place, finds place names. Type
Jagraonand pick a suggestion. - Open the SEARCH tab at the top of the left panel, next to VISUALISE. It starts with a Product name box and a note that a location filter needs an area of interest.
- Under DATA SOURCES, tick SENTINEL-2. It opens a MSI group with a cloud cover slider that starts at 100 percent, and two boxes, L1C and L2A. Tick L2A, the Level-2A product, and leave L1C unticked. Drag the slider down to 10 to 20 percent for a first try. The Filters link beside each data source opens more options, such as the satellite (S2A, S2B or S2C) and the relative orbit number.
- Define your area of interest. The pentagon button on the right of the map opens the drawing tools: upload a file, draw a rectangle, draw a polygon, or place a point. Only products that cover the area are returned.
- Scroll the left panel down to TIME RANGE. Set From and Until, and tick Filter by months if you want the same months of several years.
- Click the green Search button. Another way to get a list is to zoom the map to your area in the VISUALISE tab and click Find products for current view. The results appear under a Search link, with a line such as Showing 50 results of 503, a Select all box, an Add to workspace button and a Load more button at the bottom. Each result is a card with a thumbnail, the product name, Mission, Instrument, Size and Sensing time. The cloud percentage is not on the card. The cloud slider in step 4 already limited the results, and the thumbnail shows what the scene looks like.
- The icons on a card are an information icon, a zoom icon and a download icon. The information icon opens a Product info window with the attributes, a preview, the footprint, a Download link with product id, and the buttons Workspace and Download. Check the Size before you download: a scene at the edge of the satellite's swath can be as small as 15 MB and hold only a sliver of ground, while a full tile is about 1 GB. Click Download, or the download icon on the card, and unzip the
.zipyou receive.
The labels in steps 1 to 8 were read on the live site. The download itself was not run. A product name such as S2A_MSIL2A_20261026T101031_N0513_R022_T33UWT_20261026T150607.SAFE gives the satellite (S2A), the product level (MSIL2A), the sensing date and time, the processing baseline (N0513), the relative orbit (R022) and the tile (T33UWT). The site changes from time to time, so the labels may differ slightly from what you see.
The unzipped folder ends in .SAFE. The bands are inside GRANULE, in a folder named after the scene, then IMG_DATA, where the files are sorted by resolution:
R10mholds the 10 m bands: B02 (blue), B03 (green), B04 (red) and B08 (near infrared).R20mholds the 20 m bands, including B11 and B12 (shortwave infrared) and the scene classification layerSCL.R60mholds the 60 m bands.
The files are JPEG 2000 (.jp2). Rasterio and QGIS both open them. The values are whole numbers, not reflectance. Section "Turn the numbers into reflectance" below shows the conversion.
Route B: Python, no account
Route B asks a catalogue for the scenes, then reads only the pixels inside your area from a file that sits online. A catalogue of this kind follows a standard called STAC (SpatioTemporal Asset Catalog). Earth Search is an open STAC catalogue of Sentinel-2 data on Amazon's servers, and it needs no login.
The files there are Cloud Optimized GeoTIFFs. A Cloud Optimized GeoTIFF is a GeoTIFF organised so that a program can read one small window without downloading the rest of the file. That is why the script below downloads a few megabytes, not a gigabyte.
You need Python with requests, numpy, rasterio and shapely.
import requests
# Step 1: describe what you want. The box is west, south, east, north in longitude and latitude.
study_area = [75.40, 30.74, 75.55, 30.84] # a box of about 14 km by 11 km near Jagraon, Punjab
search = {
"collections": ["sentinel-2-c1-l2a"], # Level-2A surface reflectance
"bbox": study_area,
"datetime": "2023-02-01T00:00:00Z/2023-02-28T23:59:59Z",
"query": {"eo:cloud_cover": {"lt": 20}}, # less than 20 percent cloud over the tile
"limit": 20,
}
# Step 2: ask the catalogue. It answers with a list of scenes.
url = "https://earth-search.aws.element84.com/v1/search"
scenes = requests.post(url, json=search, timeout=60).json()["features"]
# Step 3: the catalogue returns every scene that touches your box. Keep only the scenes that cover all of it.
from shapely.geometry import shape, box
area = box(*study_area)
covering = [s for s in scenes if shape(s["geometry"]).contains(area)]
print(len(scenes), "scenes touch the box,", len(covering), "cover all of it")
# Step 4: print the date, the cloud percentage and the name of every scene that covers it
for scene in covering:
print(scene["properties"]["datetime"][:10], round(scene["properties"]["eo:cloud_cover"], 1), scene["id"])
5 scenes touch the box, 5 cover all of it
2023-02-26 0.5 S2B_T43REQ_20230226T054202_L2A
2023-02-21 6.4 S2A_T43REQ_20230221T054839_L2A
2023-02-11 0.1 S2A_T43REQ_20230211T053946_L2A
2023-02-06 9.7 S2B_T43REQ_20230206T054847_L2A
2023-02-01 0.3 S2A_T43REQ_20230201T054044_L2A
Step 3 matters. A scene is returned when its footprint touches your box, and touching is not the same as covering. Neighbouring tiles overlap by about 10 km, so a tile that only clips the edge of your box can come back with very little cloud, and a window read from it would be mostly empty. The shapely test keeps the scenes whose footprint contains the whole box, and the first line shows how many were dropped.
Each line is one scene: the date, the cloud percentage over the whole tile, and the scene name. The name tells you a lot. S2A is the satellite, T43REQ is the tile, and the number that follows is the date and time of the pass.
Pick the clearest scene and look at what the catalogue offers for it.
# Step 1: take the scene with the lowest cloud cover among those that cover your box
scene = min(covering, key=lambda s: s["properties"]["eo:cloud_cover"])
print("chosen:", scene["id"], "| cloud:", round(scene["properties"]["eo:cloud_cover"], 1), "percent")
# Step 2: every file of the scene is an asset. List their names.
print(sorted(scene["assets"].keys()))
# Step 3: look at one asset. The link points to a file you can read directly.
red = scene["assets"]["red"]
print(red["href"])
print(red["raster:bands"][0])
chosen: S2A_T43REQ_20230211T053946_L2A | cloud: 0.1 percent
['aot', 'blue', 'cloud', 'coastal', 'granule_metadata', 'green', 'nir', 'nir08', 'nir09', 'preview', 'product_metadata', 'red', 'rededge1', 'rededge2', 'rededge3', 'scl', 'snow', 'swir16', 'swir22', 'thumbnail', 'tileinfo_metadata', 'visual', 'wvp']
https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/43/R/EQ/2023/2/S2A_T43REQ_20230211T053946_L2A/B04.tif
{'nodata': 0, 'data_type': 'uint16', 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1}
The asset called red is band B04 and nir is B08. Each asset has a link and a few facts about the file. The scale and offset in raster:bands are the numbers you need for reflectance.
Turn the numbers into reflectance
Sentinel-2 stores reflectance as whole numbers to save space. To get real reflectance, which runs from 0 to about 1, you multiply by the scale and add the offset. Scenes processed with the newer baseline (04.00 and later) have an offset of -0.1, and older scenes have none. Do not remember this. Read both numbers from the catalogue each time, as the script below does.
Read only your study area
import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds
def read_window(asset, bbox):
# Step 1: open the file online. Nothing is downloaded yet.
with rasterio.open(asset["href"]) as src:
# Step 2: your box is in longitude and latitude, the file is in metres (a UTM zone), so convert the box first
left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *bbox)
window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()
# Step 3: read just that window. Only these pixels travel over the network.
values = src.read(1, window=window).astype("float32")
transform = src.window_transform(window)
crs = src.crs
# Step 4: whole numbers to reflectance. Cells with the value 0 have no data.
info = asset["raster:bands"][0]
reflectance = np.where(values == 0, np.nan, values * info["scale"] + info["offset"])
return reflectance, transform, crs
red_band, transform, crs = read_window(scene["assets"]["red"], study_area)
nir_band, _, _ = read_window(scene["assets"]["nir"], study_area)
print("window:", red_band.shape[1], "x", red_band.shape[0], "cells | CRS:", crs.to_string())
print("cells with no data:", round(float(np.isnan(red_band).mean()) * 100, 1), "percent")
print("red reflectance: from", round(float(np.nanmin(red_band)), 3), "to", round(float(np.nanmax(red_band)), 3))
print("NIR reflectance: from", round(float(np.nanmin(nir_band)), 3), "to", round(float(np.nanmax(nir_band)), 3))
window: 1440 x 1114 cells | CRS: EPSG:32643
cells with no data: 0.0 percent
red reflectance: from 0.016 to 0.454
NIR reflectance: from 0.054 to 0.645
The box is about 14 km wide and 11 km high, so the window has 1440 by 1114 cells of 10 m. Reflectance should fall between 0 and roughly 1. If you see values in the thousands, the scale was not applied. If you see negative values, check the offset.
Save what you read
# Step 1: write the two bands to one GeoTIFF with the same position and CRS as the window
with rasterio.open(
"sentinel2_red_nir.tif", "w", driver="GTiff",
height=red_band.shape[0], width=red_band.shape[1], count=2,
dtype="float32", crs=crs, transform=transform, nodata=np.nan, compress="deflate",
) as dst:
dst.write(red_band, 1)
dst.write(nir_band, 2)
dst.set_band_description(1, "red B04")
dst.set_band_description(2, "nir B08")
# Step 2: read it back to check
with rasterio.open("sentinel2_red_nir.tif") as check:
print("saved:", check.count, "bands |", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg())
saved: 2 bands | 1440 x 1114 cells | CRS EPSG 32643
Open sentinel2_red_nir.tif in QGIS and it lands on the right place, with no extra setup, because the file carries its position and CRS.
The scene classification layer
For clouds you also want the scl asset. It is a 20 m raster whose values are classes, for example 3 for cloud shadow, 8 and 9 for medium and high probability cloud, and 10 for thin cirrus. The NDVI topic shows how to use it. Read it with the same function, but remember that it has a 20 m cell and not a 10 m cell, so its window has half as many rows and columns.
Which route to choose
Route A is the better choice when you want a whole scene, want every band, or want to keep the original product with its metadata. Route B is the better choice when you want a small area, many dates, or a result you can run again next month. Route B also has no account and no quota to manage.
Common mistakes
- Downloading Level-1C because it is listed first, and then wondering why the haze never goes away.
- Judging a scene by its tile cloud percentage and skipping the preview.
- Choosing a scene that only touches your box. Keep the ones that cover all of it, and check the share of cells with no data.
- Forgetting the scale and the offset, and getting NDVI values that are wrong without any error message.
- Reading a window with longitude and latitude numbers against a file in metres.
- Mixing bands of different resolutions without resampling.
- Downloading ten full scenes, around 10 GB, for a study area that fits in a few megabytes.

