Auteur
Een stap-voor-stap handleiding voor het verzamelen en voorbewerken van satellietbeelden data voor algoritmen voor machinaal leren
Dit artikel is het eerste deel van een 2-delige serie over het gebruik van algoritmen voor machinaal leren op satellietbeelden. Hier richten we ons op de stappen data voor verzamelen en voorbewerken, terwijl we in het tweede deel het gebruik van verschillende algoritmen voor machinaal leren laten zien.
TL;DR
Dit artikel zal:
We zullen de detectie en classificatie van landbouwvelden als voorbeeld gebruiken. Python en jupyter-notebook worden gebruikt voor de voorbewerkingsstappen. Het jupyter-notebook voor het reproduceren van de stappen is beschikbaar op github.
Dit artikel gaat uit van basiskennis van data wetenschap en python.
Computer vision is een fascinerende tak van machine learning die zich de afgelopen jaren snel heeft ontwikkeld. Met openbare toegang tot veel algoritmabibliotheken(Tensorflow, OpenCv, Pytorch, Fastai...) en open source datasets(CIFAR -10, Imagenet, IMDB...) is het heel eenvoudig voor een data wetenschapper om praktijkervaring op te doen met dit onderwerp.
Maar wist je dat je ook toegang kunt krijgen tot open source satellietbeelden? Dit soort beelden kan worden gebruikt in een groot aantal toepassingen (landbouw, logistiek, energie ...). Ze bieden echter extra uitdagingen voor een data wetenschapper door hun omvang en complexiteit. Ik heb dit artikel geschreven om enkele van mijn belangrijkste ervaringen te delen.
Stap 1 - Selecteer de satelliet open data bron
Eerst moet je kiezen welke satelliet je wilt gebruiken. De belangrijkste kenmerken waar je hier op moet letten zijn :
Hier zijn satellieten waarvan data gratis beschikbaar is:
Benchmark van satellietsystemen
Aangezien het ons doel is om landbouwgewassen te detecteren en te identificeren, zou het zichtbare spectrum (dat wordt gebruikt om standaard kleurenbeelden te leveren) heel nuttig moeten zijn. Bovendien zijn er correlaties tussen het stikstofniveau van gewassen en NIR-spectrumbanden (nabij infrarood). Wat de vrije opties betreft, is Sentinel 2 momenteel de beste in termen van resolutie met 10 m voor deze spectrale banden, dus we gaan deze satelliet gebruiken.
Natuurlijk is een resolutie van 10 m veel lager dan de allernieuwste satellietbeelden van commerciële satellieten met een resolutie van 30 cm. Toegang tot commerciële satellieten data kan veel gebruikstoepassingen ontsluiten, maar ondanks de toenemende concurrentie met nieuwe spelers die deze markt betreden, blijven ze hoog geprijsd en dus moeilijk te integreren in een gebruikstoepassing op schaal.
Stap 2 - Verzamel de satellietbeelden data van Sentinel 2
Nu we onze satellietbron hebben geselecteerd, is de nieuwe stap het downloaden van de beelden. Sentinel 2 maakt deel uit van het Copernicus-programma van het Europees Ruimteagentschap. We hebben toegang tot de data met de SentinelSat python API.
We moeten eerst een account aanmaken bij de Copernicus Open Access Hub. Merk op dat deze website ook een manier biedt om afbeeldingen te bevragen met een GUI.
Een manier om de zones die je wilt verkennen te specificeren is om de API te voorzien van een geojson-bestand, wat een json is met de GPS-coördinaten positie (lengte- en breedtegraad) van een zone. Geojson kan snel worden aangemaakt op deze website.
Op deze manier kun je alle satellietbeelden opvragen die de opgegeven zone doorsnijden.
Hier geven we een lijst van alle afbeeldingen die beschikbaar zijn voor een bepaalde zone en een bepaald tijdbereik:
api = SentinelAPI(
referenties["gebruikersnaam"],
referenties["wachtwoord"],
"https://scihub.copernicus.eu/dhus"
)
shape = geojson_to_wkt(read_geojson(geojson_path))
afbeeldingen = api.query(
vorm,
date=(datum(2020, 5, 1), datum(2020, 5, 10)),
platformnaam="Sentinel-2″,
processinglevel="Level-2A",
wolkenpercentage=(0, 30)
)
images_df = api.to_dataframe(images)
Let op het gebruik van argumenten:
beelden_df
Images_df is een pandas dataframe met alle afbeeldingen die voldoen aan onze zoekopdracht. We selecteren er een met een laag cloud en downloaden deze met behulp van de api :
api.download(uuid)
Stap 3 - Begrijpen welke spectrale banden te gebruiken en kwantisatie beheren
Zodra het zipbestand is gedownload, kunnen we zien dat het eigenlijk verschillende afbeeldingen bevat.
Dit komt omdat, zoals we in stap 1 al zeiden, Sentinel 2 veel spectrale banden vastlegt en niet alleen rood, groen en blauw, en elke band in één beeld wordt vastgelegd.
Ook hebben de instrumenten niet allemaal dezelfde ruimtelijke resolutie, dus is er een submap voor elke verschillende resolutie 10, 20 en 60 meter.
If we check the GRANULE/<image_id>/IMG_DATA folder we see that there is a subfolder for each different resolution: 10, 20 and 60 meters as the spectral bands. We are going to use those with max resolution 10m as they provide Band 2, 3, 4 and 8 which are Blue, Green, Red and NIR. Band 7, 8, 9 which are vegetation red edge (spectral bands between red and NIR that mark a transition in reflectance value for vegetation) could also be interesting for detecting crops but since they are lower resolution we are going to leave them aside now
Hier laden we elke band als een 2D numpy matrix:
def get_band(image_folder, band, resolution=10):
subfolder = [f for f in os.listdir(image_folder + “/GRANULE”) if f[0] == “L”][0]
image_folder_path = f”{image_folder}/GRANULE/{subfolder}/IMG_DATA/R{resolution}m”
image_files = [im for im in os.listdir(image_folder_path) if im[-4:] == “.jp2”]
selected_file = [im for im in image_files if im.split(“_”)[2] == band][0]
with rasterio.open(f”{image_folder_path}/{selected_file}”) as infile:
img = infile.read(1)
retourneer img
band_dict = {}
voor band in ["B02", "B03", "B04", "B08"]:
band_dict[band] = get_band(image_folder, band, 10)
Een snelle test om er zeker van te zijn dat alles werkt, is door te proberen de afbeelding met echte kleuren (d.w.z. RGB) opnieuw te maken en weer te geven met opencv en matplotlib:
img = cv2.merge((band_dict["B04"], band_dict["B03"], band_dict["B02"]))
plt.imshow(img)
Satellietbeeld - kwestie van kleurintensiteit (Copernicus Sentinel data 2020)
Er is echter een probleem met de kleurintensiteit. Dit komt doordat de kwantisatie (aantal mogelijke waarden) afwijkt van de standaard voor matplotlib die uitgaat van een 0-255 int waarde of een 0-1 float waarde. Standaardafbeeldingen gebruiken vaak een 8-bits kwantisatie (256 mogelijke waarden), maar Sentinel 2-afbeeldingen gebruiken een 12-bits kwantisatie en een herbewerking converteert de waarden naar een 16-bits geheel getal (65536 waarden).
Als we een 16-bits geheel getal willen schalen naar een 8-bits geheel getal, moeten we in theorie delen door 256, maar dit resulteert in een erg donkere afbeelding omdat niet het hele bereik van mogelijke waarden wordt gebruikt. De maximale waarde in onze afbeelding is eigenlijk 16752 van de 65536 en weinig pixels bereiken waarden hoger dan 4000. Delen door 8 geeft dus een afbeelding met een behoorlijk contrast.
img_verwerkt = img / 8
img_bewerkt = img_bewerkt.astype(int)
plt.imshow(img_bewerkt)
Satellietbeeld - Kleurintensiteit rechts (Copernicus Sentinel data 2020)
We hebben nu een 3-dimensionale tensor van de vorm (10980, 10980, 3) van 8-bits gehele getallen in de vorm van een numpy array
Merk op dat :
Stap 4 - Splits de afbeeldingen op naar machine learning-klare formaten
We zouden kunnen proberen een algoritme toe te passen op onze huidige afbeelding, maar in de praktijk zou dat met de meeste technieken niet werken vanwege de grootte. Zelfs met veel rekenkracht zouden de meeste algoritmen (en vooral deep learning) niet werken.
We moeten de afbeelding dus opsplitsen in fragmenten. Een vraag is hoe we de grootte van het fragment kunnen instellen.
We besluiten om elke afbeelding op te delen in een raster van 45 * 45 (45 is gemakshalve een deler van 10980)
Voorbij deze parameterbepaling is deze stap rechttoe rechtaan met behulp van numpy array manipulatie . We slaan onze waarden op in een dict met (x, y) tupels als sleutels.
frag_count = 45
frag_size = int(img_processed.shape[0] / frag_count)
frag_dict = {}
voor y, x in itertools.product(range(frag_count), range(frag_count)):
frag_dict[(x, y)] = img_processed[y*frag_size: (y+1)*frag_size,
x*frag_size: (x+1)*frag_size, :]
plt.imshow(frag_dict[(10, 10)])
Fragment van het satellietbeeld (Copernicus Sentinel data 2020)
Stap 5 - Koppel je data aan de satellietbeelden met behulp van GPS-naar-UTM-conversie
Tot slot moeten we onze afbeeldingen koppelen aan de data die we hebben over landbouwvelden om een gesuperviseerde machine-learning aanpak mogelijk te maken. In ons geval hebben we data over landbouwvelden in de vorm van een lijst GPS-coördinaten (breedtegraad, lengtegraad), vergelijkbaar met de geojson die we in stap 2 hebben gebruikt en we willen hun pixelpositie kunnen vinden. Voor de duidelijkheid zal ik echter eerst een pixelconversie naar GPS-coördinaten laten zien en daarna het omgekeerde
De metadata die zijn verzameld tijdens het opvragen van de SentinelSat API geeft de GPS-coördinaten van de hoeken van de afbeelding (footprint kolom) en we willen de coördinaten van een specifieke pixel. Breedtegraad en lengtegraad, die hoekwaarden zijn, evolueren niet lineair binnen een afbeelding, dus we kunnen geen eenvoudige verhouding maken met behulp van de pixelpositie. Er zijn manieren om dit op te lossen met wiskundige vergelijkingen, maar er bestaan al kant-en-klare oplossingen in python.
Een snelle oplossing is om de pixelpositie te converteren naar UTM (Universal Transverse Mercator), waarin een positie wordt gedefinieerd door een zone en een (x, y) positie (in meter-eenheid) die lineair evolueert met de afbeelding. We kunnen dan een UTM-naar- breedtegraad/lengtegraad conversie gebruiken in de utm bibliotheek. Om dit te doen moeten we eerst de UTM-zone en positie van de linkerbovenhoek van de afbeelding opvragen. Deze kunnen worden gevonden in de metagegevens van de afbeelding.
#Vleesgegevens verkrijgen
transform = rasterio.open(tci_file_path, driver='JP2OpenJPEG').transform
zone_nummer = int(tci_file_path.split("/")[-1][1:3])
zone_letter = tci_file_path.split("/")[-1][0]
utm_x, utm_y = transform[2], transform[5]
# Pixelpositie omzetten naar utm
oost = utm_x + pixel_kolom * 10
noord = utm_y + pixel_rij * - 10
# UTM omzetten naar breedtegraad en lengtegraad
breedtegraad, lengtegraad = utm.to_latlon(east, north, zone_nummer, zone_letter)
De omzetting van GPS-positie naar pixel gebeurt met de omgekeerde formules. In dat geval moeten we ervoor zorgen dat de verkregen zone overeenkomt met onze afbeelding. Dit kan worden gedaan door de zone te specificeren. Als ons object niet aanwezig is in de afbeelding, resulteert dit in een pixelwaarde die buiten de grenzen valt.
# Breedtegraad en lengtegraad omzetten naar UTM
oost, noord, zone_nummer, zone_letter = utm.from_latlon(
breedtegraad, lengtegraad, force_zone_number=zone_number
)
# UTM omzetten naar kolom en rij
pixel_kolom = rond((oost - utm_x) / 10)
pixel_row = rond((noord - utm_y) / -10)
onbewerkt bekijken
Conclusie
We zijn nu klaar om de satellietbeelden te gebruiken voor machine learning!
In het volgende artikel van deze 2-delige serie zullen we zien hoe we die beelden kunnen gebruiken om veldoppervlakken te detecteren en classificeren met behulp van zowel supervised als unsupervised machine learning.
Bedankt voor het lezen en aarzel niet om de Artefact tech blog te volgen als je op de hoogte wilt worden gebracht wanneer dit volgende artikel verschijnt!
