diff --git a/Acolite-atmospheric-correction.ipynb b/Acolite-atmospheric-correction.ipynb new file mode 100644 index 0000000..6c7fc73 --- /dev/null +++ b/Acolite-atmospheric-correction.ipynb @@ -0,0 +1,459 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import glob\n", + "import time\n", + "import os\n", + "import re\n", + "from osgeo import gdal\n", + "import warnings\n", + "from pyproj import Proj, transform\n", + "import osr\n", + "import numpy as np\n", + "import shutil as sh\n", + "from affine import Affine\n", + "import shutil\n", + "import geopandas\n", + "import rasterio\n", + "import utm\n", + "from shapely.geometry import Point\n", + "warnings.filterwarnings(\"ignore\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def save_multiband (output_name, dataset, raster_data, driver,\n", + " NaN_Value, nband, option, arr_projection=None):\n", + " if arr_projection is None:\n", + " arr_projection = []\n", + " else:\n", + " if str(type(arr_projection)) == \"\":\n", + " srs = arr_projection.GetProjectionRef()\n", + " elif str(type(arr_projection)) == \"\":\n", + " srs = osr.SpatialReference()\n", + " srs.ImportFromEPSG(arr_projection)\n", + " srs = srs.ExportToPrettyWkt()\n", + " else:\n", + " srs = arr_projection\n", + " NaN_rast = NaN_Value\n", + " # Open the reference dataset\n", + " g = dataset\n", + " # Get the Geotransform vector\n", + " if raster_data is False:\n", + " raster_data = g.ReadAsArray()\n", + " rast_arr = np.copy(raster_data)\n", + " #auxiliar\n", + " band = 0\n", + " if type(raster_data) == tuple:\n", + " rast_arr = np.array(raster_data[band])\n", + " if str(type(g)) == \"\":\n", + " geo_transform = g.GetGeoTransform()\n", + " x_size = g.RasterXSize # Raster xsize\n", + " y_size = g.RasterYSize # Raster ysize\n", + " srs = g.GetProjectionRef() # Projection\n", + " elif str(type(g)) == \"\":\n", + " geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])\n", + " rast_arr = raster_data[band,:,:]\n", + " x_size = int(rast_arr.shape[1])\n", + " y_size = int(rast_arr.shape[0])\n", + " \n", + " g = dataset\n", + " driver = gdal.GetDriverByName(\"GTiff\")\n", + " dataset_out = driver.Create(output_name, x_size, y_size, nband,\n", + " gdal.GDT_Float32)\n", + " #rast_arr[rast_arr == NaN_rast] = np.NaN\n", + " geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])\n", + " #dataset_out.SetGeoTransform(geo_transform)\n", + " #dataset_out.SetProjection(srs)\n", + " #end auxiliar\n", + " for band in range(nband):\n", + " if type(raster_data) == tuple:\n", + " rast_arr = np.array(raster_data[band])\n", + " if str(type(g)) == \"\":\n", + " geo_transform = g.GetGeoTransform()\n", + " x_size = g.RasterXSize # Raster xsize\n", + " y_size = g.RasterYSize # Raster ysize\n", + " srs = g.GetProjectionRef() # Projection\n", + " elif str(type(g)) == \"\":\n", + " geo_transform = (g[2], g[0], g[1], g[5], g[3], g[4])\n", + " rast_arr = raster_data[band,:,:]\n", + " x_size = int(rast_arr.shape[1])\n", + " y_size = int(rast_arr.shape[0])\n", + " #PROCESS RASTERIO NUMPY\n", + " else:\n", + " geo_transform = (g[1][2], g[1][0], g[1][1], g[1][5], g[1][3], g[1][4])\n", + " rast_arr = np.array(g[0])\n", + " x_size = int(rast_arr.shape[2])\n", + " y_size = int(rast_arr.shape[1])\n", + "\n", + " if option==1:\n", + " rast_arr[rast_arr == NaN_rast] = np.NaN\n", + " dataset_out.SetGeoTransform(geo_transform)\n", + " dataset_out.SetProjection(srs)\n", + " dataset_out.GetRasterBand(band + 1).WriteArray(\n", + " rast_arr.astype(np.float32))\n", + " elif option==2:\n", + " dataset_out.SetGeoTransform(geo_transform)\n", + " dataset_out.SetProjection(srs)\n", + " dataset_out.GetRasterBand(band + 1).WriteArray(\n", + " rast_arr.astype(np.float32))\n", + " \n", + "# rast_arr[rast_arr == NaN_rast] = np.NaN \n", + "# dataset_out.SetGeoTransform(geo_transform)\n", + "# dataset_out.SetProjection(srs)\n", + "# dataset_out.GetRasterBand(band + 1).WriteArray(rast_arr.astype(np.float32))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def GetTransform(Transform):\n", + " transform = Transform.GetGeoTransform()\n", + " Aff = Affine(transform[1], transform[2], transform[0], transform[4],\n", + " transform[5], transform[3])\n", + " return(Aff)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "#Get directory with all raw sat images for Rayleigh atmospheric correction\n", + "path_to_s2 = r'/mnt/data4/antonis/dataset_new'\n", + "s2_files = glob.glob(path_to_s2 + \"/*\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "#Step1: Specify product name to be ACOLITE corrected" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801 is present in the list\n" + ] + } + ], + "source": [ + "requested_product = 'S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801'\n", + "for k in range(len(s2_files)):\n", + " if requested_product in s2_files[k][-60:]:\n", + " print('%s is present in the list' % requested_product)\n", + " requested_path = s2_files[0][:-60] + requested_product" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'/mnt/data4/antonis/dataset_new/S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801/RAW/S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801.SAFE'" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#Step2: Show path of RAW data folder to ACOLITE correction\n", + "requested_path+'/RAW/'+requested_product+'.SAFE'" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Step1: Copying settings_file from acolite to image folder\n", + "Step2: Creation of folder of the desired atmosphericly corrected product\n", + "Successfully created the directory /mnt/data4/antonis/acolite_corrected/S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801.SAFE \n", + "Now processing: S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801\n", + "Step3: Change inputs inside settings file (inputfile, outputfile, extend)\n", + "/mnt/data4/antonis/dataset_new/S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801/RAW/S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801.SAFE/GRANULE/L1C_T39RUM_A007039_20180712T072748/IMG_DATA/T39RUM_20180712T071619_B01.jp2\n", + "Step4: Update extend parameter in settings_file using gdal\n", + "Step5: Run ACOLITE based on the updated settings_file\n", + "Running... /home/antonis/acolite_py_linux/dist/acolite/acolite --cli --settings=\"/home/antonis/sentinel-2/settings_file.txt\"\n", + "Acolite algorithm completed for: S2B_MSIL1C_20180712T071619_N0206_R006_T39RUM_20180712T110801 in: 0:31:26.04\n" + ] + } + ], + "source": [ + "#Step3: Run ACOLITE\n", + "#Step 1:\n", + "print('Step1: Copying settings_file from acolite to image folder')\n", + "#Copy settings_file from acolite folder to image folder to proceed atmospherical corrections\n", + "src = r'/home/antonis/acolite_py_linux/dist/acolite/settings_file.txt'\n", + "dst =r'/home/antonis/sentinel-2/'\n", + "sh.copy(src, dst)\n", + "\n", + "#Step2:\n", + "print('Step2: Creation of folder of the desired atmosphericly corrected product')\n", + "#Create folder of the desired atmosphericly corrected product\n", + "output_folder = r'/mnt/data4/antonis/acolite_corrected' +str('/') + requested_product +'.SAFE'\n", + "if not os.path.exists(output_folder):\n", + " os.mkdir(output_folder)\n", + " print (\"Successfully created the directory %s \" % output_folder)\n", + "else:\n", + " print (\"Directory %s already exists\" % output_folder)\n", + "\n", + "print('Now processing: ',requested_product)\n", + "\n", + "#Step3:\n", + "print('Step3: Change inputs inside settings file (inputfile, outputfile, extend)')\n", + "#Change inputs inside settings file (inputfile, outputfile, extend)\n", + "#Currently parameter specific based on the original settings_file located in the acolite folder\n", + "file = r'/home/antonis/sentinel-2/settings_file.txt'\n", + "ouput = r'/home/antonis/sentinel-2/settings_file_n.txt'\n", + "output_dir = output_folder\n", + "\n", + "with open(file, 'r') as f:\n", + " with open(ouput, 'wt') as new_f:\n", + " for line in f:\n", + " new_line = line.replace(r'inputfile=C:\\Users\\antonkout\\Documents\\Work\\NTUA\\Data\\Preprocess\\LE07_L2SP_021040_20100501_20200911_02_T1', 'inputfile='+requested_path+'/RAW/'+requested_product+'.SAFE')\n", + " new_f.write(new_line)\n", + "\n", + "file = r'/home/antonis/sentinel-2/settings_file_n.txt'\n", + "ouput = r'/home/antonis/sentinel-2/settings_file.txt'\n", + "\n", + "with open(file, 'r') as f:\n", + " with open(ouput, 'wt') as new_f:\n", + " for line in f:\n", + " new_line = line.replace(r'output=C:\\Users\\antonkout\\Documents\\Work\\NTUA\\Data\\Preprocess\\LE07_L2SP_021040_20100501_20200911_02_T1_results', 'output='+output_dir)\n", + " new_f.write(new_line)\n", + "\n", + "if os.path.exists(file): \n", + " os.remove(file)\n", + "else: \n", + " print(\"The file does not exist\") \n", + "\n", + "#To find extend must first find the desired sat img\n", + "for root, dirs, files in os.walk(requested_path+'/RAW/'+requested_product+'.SAFE'):\n", + " for name in files:\n", + " if name.endswith(\"_B01.jp2\"):\n", + " sat = root+str(r'/')+name\n", + " print(sat)\n", + "#Step4\n", + "print('Step4: Update extend parameter in settings_file using gdal')\n", + "#Get extend using gdal\n", + "ds = gdal.Open(sat, 1)\n", + "geotransform = ds.GetGeoTransform()\n", + "proj = osr.SpatialReference(wkt=ds.GetProjection())\n", + "inprj = 'epsg:'+proj.GetAttrValue('AUTHORITY',1)\n", + "\n", + "upx, xres, xskew, upy, yskew, yres = geotransform\n", + "cols = ds.RasterXSize\n", + "rows = ds.RasterYSize\n", + "\n", + "ulx = upx + 0*xres + 0*xskew\n", + "uly = upy + 0*yskew + 0*yres\n", + "\n", + "llx = upx + 0*xres + rows*xskew\n", + "lly = upy + 0*yskew + rows*yres\n", + "\n", + "lrx = upx + cols*xres + rows*xskew\n", + "lry = upy + cols*yskew + rows*yres\n", + "\n", + "urx = upx + cols*xres + 0*xskew\n", + "ury = upy + cols*yskew + 0*yres\n", + "\n", + "inProj = Proj(init=inprj)\n", + "outProj = Proj(init='epsg:4326')\n", + "\n", + "#Transform coordinated to desired epsg\n", + "ext = (ulx,uly),(llx,lly),(lrx,lry),(urx,ury)\n", + "extend=[]\n", + "for k in range(len(ext)):\n", + " extend.append(transform(inProj,outProj,ext[k][0],ext[k][1]))\n", + "xmax,xmin,ymax,ymin = np.asarray(extend)[:,1].max(),np.asarray(extend)[:,1].min(),np.asarray(extend)[:,0].max(),np.asarray(extend)[:,0].min()\n", + "extend = xmin,ymin,xmax,ymax\n", + "extend = str(extend).replace('(','').replace(')','')\n", + "\n", + "#Replace with the new extend\n", + "file =r'/home/antonis/sentinel-2/settings_file.txt'\n", + "ouput = r'/home/antonis/sentinel-2/settings_file_n.txt'\n", + "\n", + "with open(file, 'r') as f:\n", + " with open(ouput, 'wt') as new_f:\n", + " for line in f:\n", + " new_line = line.replace(r'limit= 27.883946, -90.264583, 29.858269, -87.733565', 'limit= '+extend)\n", + " new_f.write(new_line)\n", + "\n", + "#Rename text file to the origin name and remove the temp one\n", + "if os.path.exists(file): \n", + " os.remove(file)\n", + "else: \n", + " print(\"The file does not exist\") \n", + "os.rename(ouput,file)\n", + "\n", + "#Step 5:\n", + "print('Step5: Run ACOLITE based on the updated settings_file')\n", + "#Run ACOLITE based on the constructed settings_file\n", + "starttime=time.time()\n", + "cmd = r'/home/antonis/acolite_py_linux/dist/acolite/acolite --cli --settings='+'\"{}\"'.format(file)\n", + "print (\"Running...\", cmd)\n", + "os.system(cmd)\n", + "elapsedtime=time.time()-starttime\n", + "mins, secs = divmod(elapsedtime, 60)\n", + "hours, mins = divmod(mins, 60)\n", + "print (\"Acolite algorithm completed for:\", requested_product, \"in: {}:{}:{}\".format(int(hours), int(mins), round(secs,2)))\n", + "\n", + "#Remove settings_file in order to copy it again from the acolite installed folder\n", + "#os.remove(file)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "#Step4: Create multiband rhorc\n", + "sat = []\n", + "indx = ['_443.tif', '_492.tif', '_560.tif', '_665.tif', '_704.tif', '_740.tif', '_783.tif', '_833.tif', '_865.tif', '_1614.tif', '_2202.tif']\n", + "indxn = ['_442.tif', '_492.tif', '_559.tif', '_665.tif', '_704.tif', '_739.tif', '_780.tif', '_833.tif', '_864.tif', '_1610.tif', '_2186.tif']\n", + "\n", + "#To find extend must first find the desired sat img\n", + "for root, dirs, files in os.walk(output_folder):\n", + " for name in files:\n", + " if 'rhorc' in name:\n", + " sat.append(root+str(r'/')+name)\n", + "rhorc_sat=[]\n", + "for a in range(len(sat)):\n", + " if any(sat[a][-8:] in s for s in indx):\n", + " flag = 1\n", + " rhorc_sat.append(sat[a])\n", + "\n", + "if len(rhorc_sat)<10:\n", + " rhorc_sat = []\n", + " flag = 2 \n", + " for a in range(len(sat)):\n", + " if any(sat[a][-8:] in s for s in indxn):\n", + " rhorc_sat.append(sat[a])\n", + "\n", + "rhorc_index =[]\n", + "for k in range(len(rhorc_sat)):\n", + " if flag == 1:\n", + " rhorc_index.append(rhorc_sat[k][:143]+indx[k])\n", + " elif flag == 2:\n", + " rhorc_index.append(rhorc_sat[k][:143]+indxn[k])\n", + " \n", + "B1, B2, B3, B4, B5, B6, B7, B8, B9, B10, B11 = gdal.Open(rhorc_index[0]), gdal.Open(rhorc_index[1]), gdal.Open(rhorc_index[2]), gdal.Open(rhorc_index[3]),gdal.Open(rhorc_index[4]), gdal.Open(rhorc_index[5]), gdal.Open(rhorc_index[6]), gdal.Open(rhorc_index[7]), gdal.Open(rhorc_index[8]), gdal.Open(rhorc_index[9]), gdal.Open(rhorc_index[10])\n", + "srs = B1.GetProjectionRef()\n", + "b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11 = (B1.ReadAsArray(), B2.ReadAsArray(),B3.ReadAsArray(), B4.ReadAsArray(),\n", + " B5.ReadAsArray(), B6.ReadAsArray(),B7.ReadAsArray(), B8.ReadAsArray(),\n", + " B9.ReadAsArray(), B10.ReadAsArray(), B11.ReadAsArray())\n", + "B = np.stack([b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11])\n", + "for k in range(B.shape[0]):\n", + " B[k] = np.where(B[k]==9.96921e+36, 0, B[k])\n", + " \n", + "rhorc_file = requested_path+'/ACOLITE/'+requested_product+'_rhorc.tif'\n", + "option = 1\n", + "save_multiband(rhorc_file, GetTransform(B1), B, \"GTiff\", 0, B.shape[0],option, srs)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "#Step5: Delete folder of acolite corrected\n", + "shutil.rmtree(output_folder)" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "#Step6: Create RGB composite\n", + "sat = []\n", + "indx = ['_B02.jp2', '_B03.jp2', '_B04.jp2']\n", + "#To find extend must first find the desired sat img\n", + "for root, dirs, files in os.walk(requested_path+'/RAW/'+requested_product+'.SAFE'):\n", + " for name in files:\n", + " if '.jp2' in name:\n", + " sat.append(root+str(r'/')+name)\n", + "raw_sat=[]\n", + "for a in range(len(sat)):\n", + " if any(sat[a][-8:] in s for s in indx):\n", + " raw_sat.append(sat[a])\n", + "\n", + "rgb_index =[]\n", + "for k in range(len(raw_sat)):\n", + " rgb_index.append(raw_sat[k][:-8]+indx[k])\n", + " \n", + "B2, B3, B4 = gdal.Open(rgb_index[0]), gdal.Open(rgb_index[1]), gdal.Open(rgb_index[2])\n", + "srs = B2.GetProjectionRef()\n", + "b2, b3, b4 = (B2.ReadAsArray(), B3.ReadAsArray(),B4.ReadAsArray())\n", + "B = np.stack([b2, b3, b4])\n", + "\n", + "rgb_file = requested_path+'/RGB/'+requested_product+'_rgb.tif'\n", + "option = 2\n", + "save_multiband(rgb_file, GetTransform(B2), B, \"GTiff\", 0, B.shape[0], option, srs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}