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 zullen we ons richten op de stappen van de data verzameling en voorbewerking, terwijl het tweede deel het gebruik van verschillende algoritmen voor machinaal leren zal 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. De 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 machinaal leren die zich de afgelopen jaren snel heeft ontwikkeld. Met openbare toegang tot veel bibliotheken met algoritmen (Tensorflow, OpenCv, Pytorch, Fastai ...) evenals open source datasets (CIFAR -10, Imagenet, IMDB ...) is het heel gemakkelijk voor een data wetenschapper om praktijkervaring met dit onderwerp op te doen.
Maar wist u dat u ook toegang kunt krijgen tot open source satellietbeelden? Dergelijke beelden kunnen in een groot aantal toepassingen gebruikt worden (Landbouw, Logistiek, Energie ...). Ze bieden echter extra uitdagingen voor een data wetenschapper vanwege 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 u kiezen welke satelliet u wilt gebruiken. De belangrijkste kenmerken waar u hier op moet letten zijn :
Hier zijn satellieten waarvan de data gratis beschikbaar zijn:
Benchmark van satellietsystemen
Aangezien het ons doel is om landbouwgewassen te detecteren en te identificeren, zou het zichtbare spectrum (dat gebruikt wordt 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 gratis opties betreft, is Sentinel 2 momenteel de beste wat resolutie betreft met 10 m voor die spectrale banden, dus gaan we deze satelliet gebruiken.
Natuurlijk is een resolutie van 10 m veel lager dan de allernieuwste satellietbeelden met commerciële satellieten die een resolutie van 30 cm bereiken. Toegang tot commerciële data-satellieten kan veel gebruikssituaties ontsluiten, maar ondanks de toenemende concurrentie met nieuwe spelers die deze markt betreden, blijven ze hoog geprijsd en dus moeilijk te integreren in een gebruikssituatie 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 de Europees Ruimteagentschap. We hebben toegang tot de data met de SentinelSat python API.
We moeten eerst een account aanmaken bij de Copernicus Open Toegang Hub. Merk ook op dat deze website een manier biedt om afbeeldingen met een GUI te bevragen.
Een manier om de zones die u wilt verkennen op te geven, is door de API te voorzien van een geojson-bestand, wat een json is die de positie van de GPS-coördinaten (breedtegraad en lengtegraad) van een zone bevat. Geojson kan hier snel worden aangemaakt website.
Op deze manier kunt u alle satellietbeelden opvragen die de opgegeven zone doorsnijden.
Hier geven we een lijst van alle afbeeldingen die beschikbaar zijn op een voor een bepaalde zone en tijdsbereik:
api = SentinelAPI(
referenties[“gebruikersnaam”],
referenties[“wachtwoord”],
“https://scihub.copernicus.eu/dhus”
)
shape = geojson_to_wkt(read_geojson(geojson_path))
afbeeldingen = api.query(
vorm,
date=(date(2020, 5, 1), date(2020, 5, 10)),
platformname=”Sentinel-2″,
verwerkingsniveau = “Niveau-2A”,
cloudbedekkingspercentage=(0, 30)
)
images_df = api.to_dataframe(images)
Let op het gebruik van argumenten:
afbeeldingen_df
Images_df is een pandas dataframe dat alle afbeeldingen bevat die overeenkomen met onze zoekopdracht. We selecteren er een met een lage cloud dekking en downloaden deze met behulp van de api :
api.download(uuid)
Stap 3 - Begrijpen welke spectrale banden u moet gebruiken en kwantisatie beheren
Zodra het zipbestand is gedownload, kunnen we zien dat het eigenlijk meerdere 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.
Als we de GRANULE//IMG_DATA map zien we dat er een submap is voor elke verschillende resolutie: 10, 20 en 60 meter als de spectrale banden. We gaan die met de maximale resolutie van 10 m gebruiken, omdat deze band 2, 3, 4 en 8 bieden, die blauw, groen, rood en NIR zijn. Band 7, 8 en 9, die de rode rand van de vegetatie vormen (spectrale banden tussen rood en NIR die een overgang in reflectiewaarde voor vegetatie markeren), zouden ook interessant kunnen zijn voor het detecteren van gewassen, maar omdat ze een lagere resolutie hebben, laten we deze nu buiten beschouwing.
Hier laden we elke band als een 2D numpy matrix:
def get_band(image_folder, band, resolutie=10):
subfolder = [f for f in os.listdir(image_folder + "/GRANULE") if f[0] == "L"][0]
image_folder_path = f"/GRANULE//IMG_DATA/Rm"
image_files = [im for im in os.listdir(image_folder_path) if im[-4:] == ".jp2"]
geselecteerd_bestand = [im for im in image_files if im.split("_")[2] == band][0]
met rasterio.open(f"/") als 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 proberen om de afbeelding in ware kleur (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 - Probleem met 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 maximumwaarde in onze afbeelding is in feite 16752 van de 65536 en er zijn maar weinig pixels die waarden hoger dan 4000 bereiken, dus delen door 8 geeft een afbeelding met een behoorlijk contrast.
img_verwerkt = img / 8
img_verwerkt = img_verwerkt.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 tot machine-learningklare formaten
We zouden kunnen proberen om een algoritme op onze huidige afbeelding toe te passen, 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 in fragmenten opsplitsen. Een vraag is hoe we de grootte van het fragment kunnen instellen.
We besluiten om elke afbeelding op te splitsen 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 arraymanipulatie . 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 uw 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-learningbenadering mogelijk te maken. In ons geval hebben we data op landbouwvelden in de vorm van een lijst met 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 pixelomzetting naar GPS-coördinaten laten zien en daarna het omgekeerde
De metadata die verzameld is tijdens het opvragen van de SentinelSat API geeft de GPS-coördinaten van de hoeken van de afbeelding (footprint kolom) en wij 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 door de pixelpositie te gebruiken. 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 met de afbeelding evolueert. We kunnen dan een UTM-naar- breedtegraad/lengtegraad conversie gebruiken die in de utm-bibliotheek beschikbaar is. Om dit te doen, moeten we eerst de UTM-zone en positie van de linkerbovenhoek van de afbeelding verkrijgen. Deze kunnen gevonden worden in de metadata van de afbeelding in ware kleur.
#Obevattend meatadata
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
noorden = utm_y + pixel_row * - 10
# UTM omrekenen naar breedtegraad en lengtegraad
breedtegraad, lengtegraad = utm.to_latlon(east, north, zone_nummer, zone_letter)
De omzetting van GPS-positie naar pixel wordt gedaan met behulp van de omgekeerde formules. In dat geval moeten we ervoor zorgen dat de verkregen zone overeenkomt met onze afbeelding. Dit kan gedaan worden door de zone op te geven. Als ons object niet in de afbeelding aanwezig is, zal dit resulteren in een pixelwaarde die buiten de grenzen valt.
# Breedtegraad en lengtegraad omrekenen naar UTM
oost, noord, zone_nummer, zone_letter = utm.from_latlon(
breedtegraad, lengtegraad, force_zone_number=zone_number
)
# UTM omzetten naar kolom en rij
pixe_column = rond((oost - utm_x) / 10)
pixel_row = rond((noord - utm_y) / -10)
rauw bekijken
Conclusie
We zijn nu klaar om de satellietbeelden te gebruiken voor machinaal leren!
In het volgende artikel van deze 2-delige serie zullen we zien hoe we die afbeeldingen 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 volg de Artefact tech blog als u op de hoogte wilt worden gebracht wanneer dit volgende artikel uitkomt!
