add file

consts.py

@@ -39,3 +39,29 @@
 OUT_GRID_POLYGON_FILE = ''
 POINT_OFFSET_TO_GENERATE_GRID = 500
 POINT_FILE_COORDINATE_SYS = 'GCS'
+
+
+# tool: clip raster by extent
+RASTER_INPUT_FILE = ''
+CLIP_BOUNDARY_FILE = ''
+CLIP_RASTER_OUTPUT_FILE = ''
+
+# tool: convert raster pixel to points
+RASTER_INPUT_FILE = ''
+OUT_POINT_FILE = ''
+
+
+# tool: convert dn to radiance of satellite image bands
+RASTER_INPUT_FILE = ''
+RASTER_WORKING_DIR = ''
+BANDS_LMAX_VALUES = [193.0, 365.0, 264.0, 221.0, 30.2, 16.5]
+BANDS_LMIN_VALUES = [-1.520, -2.840, -1.170, -1.510, -0.370, -0.150]
+QCAL_MIN_VALUE = 1
+QCAL_MAX_VALUE = 255
+OUT_NAME_PREFIX = 'band_'
+
+
+# tool: export lulc individual class
+LULC_RASTER_FILE = ''
+LULC_OUT_CLASS_FILE = ''
+LULC_CLASS_CODE = 1

docs/functions.md

@@ -12,6 +12,8 @@ flip_line
 interval_points_straightline
 offset_random_point
 read_raster_as_array
+np2gdal_dtype
+write_geotiff_file
 ```
 
 ### ref_scripts.py

src/raster_ops/components/clip_raster.py

@@ -12,17 +12,19 @@
 from osgeo import gdal, gdalnumeric, ogr
 from exception import CustomException
 from ref_scripts import world2Pixel
+from raster_ops.config_entity import RasterdataConfig
 
 
 class ClipRaster:
     """
     Clip Raster Data.
     """
 
-    def __init__(self, inRaster, inPolygon):
+    def __init__(self, rstconfig: RasterdataConfig):
         """constructure"""
-        self.inRaster = inRaster
-        self.inPolygon = inPolygon
+        self.rstconfig = rstconfig
+        self.inRaster = self.rstconfig.raster_infile
+        self.inPolygon = self.rstconfig.polygon_infile
 
     def subset_by_extent(self):
         # Open the image as a read only image

src/raster_ops/components/convert_raster_to_point.py

@@ -8,6 +8,7 @@
 """
 
 from osgeo import gdal, gdalconst, gdalnumeric, ogr, osr
+from exception import CustomException
 
 
 def raster_to_point(raster_data, outpoint_file, outfield="DN"):
@@ -18,8 +19,8 @@ def raster_to_point(raster_data, outpoint_file, outfield="DN"):
     geop = ds.GetProjection()
     min_x = geot[0]
     max_y = geot[3]
-    max_x = min_x + geot[1]*ds.RasterXSize
-    min_y = max_y + geot[5]*ds.RasterYSize
+    # max_x = min_x + geot[1]*ds.RasterXSize
+    # min_y = max_y + geot[5]*ds.RasterYSize
     # get spatial reference from raster image
     sr = osr.SpatialReference()
     sr.ImportFromWkt(ds.GetProjection())
@@ -34,25 +35,29 @@ def raster_to_point(raster_data, outpoint_file, outfield="DN"):
     lyrDef = lyr.GetLayerDefn()
     # dataset to array
     array = gdalnumeric.LoadFile(raster_data)
-    # iterate over the numpy points
-    for i in range(array.shape[0]):
-        for j in range(array.shape[1]):
-            try:
-                x = j*geot[1] + min_x + geot[1]/2.0
-                y = i*geot[5] + max_y + geot[5]/2.0
-                # create points
-                point = ogr.Geometry(ogr.wkbPoint)
-                point.SetPoint(0, x, y)
-                # put point as a geometry inside a feature
-                featIndex = 0
-                feat = ogr.Feature(lyrDef)
-                feat.SetField(outfield, array[i][j])
-                feat.SetGeometry(point)
-                feat.SetFID(featIndex)
-                # put feature in a layer
-                lyr.CreateFeature(feat)
-            except:
-                pass
+    # check array dimension
+    if array.ndim == 3:
+        raise CustomException("found multi-channel raster data")
+    else:
+        # iterate over the numpy points
+        for i in range(array.shape[0]):
+            for j in range(array.shape[1]):
+                try:
+                    x = j*geot[1] + min_x + geot[1]/2.0
+                    y = i*geot[5] + max_y + geot[5]/2.0
+                    # create points
+                    point = ogr.Geometry(ogr.wkbPoint)
+                    point.SetPoint(0, x, y)
+                    # put point as a geometry inside a feature
+                    featIndex = 0
+                    feat = ogr.Feature(lyrDef)
+                    feat.SetField(outfield, array[i][j])
+                    feat.SetGeometry(point)
+                    feat.SetFID(featIndex)
+                    # put feature in a layer
+                    lyr.CreateFeature(feat)
+                except:
+                    pass
     # Flush
     shapeData.Destroy()
-    return "Process Completed!"
+    return

src/raster_ops/components/dn_to_radiance.py

@@ -16,38 +16,30 @@ def dn_to_radiance(input_image, output_folder, lmax_list, lmin_list, QCALMIN, QC
     Digital Number to Radiance Conversion.
     """
     # set output file projection from input file
-    fopen = gdal.Open(input_image)
-    geotransform = fopen.GetGeoTransform()
-    prj = fopen.GetProjection()
+    ds = gdal.Open(input_image)
+    geotransform = ds.GetGeoTransform()
+    prj = ds.GetProjection()
 
     # get total input raster bands (count)
-    inrasband = [band+1 for band in range(fopen.RasterCount)]
+    inrasband = [band+1 for band in range(ds.RasterCount)]
 
     # set output file driver
     # HFA for *.img reading and writting
     outdriver = gdal.GetDriverByName('GTiff')
 
     # processing
     for i in inrasband:
-        dnfile = fopen.GetRasterBand(i)
+        dnfile = ds.GetRasterBand(i)
         dn = numpy.array(dnfile.ReadAsArray())
         radiance = ((lmax_list[i-1] - lmin_list[i-1]) /
                     (QCALMAX - QCALMIN))*dn + lmin_list[i-1]
         # create output file with projection
         radresult = gdal_array.OpenArray(radiance)
-        geo = radresult.SetGeoTransform(geotransform)
-        proj = radresult.SetProjection(prj)
+        radresult.SetGeoTransform(geotransform)
+        radresult.SetProjection(prj)
         # print ("Digital number to radiance conversion for band : %d" % i)
-        outputs = outdriver.CreateCopy(
-            output_folder+"\\"+prefix+"_"+str(i)+".img", radresult)
+        outdriver.CreateCopy(output_folder+"\\"+prefix +
+                             "_"+str(i)+".img", radresult)
+    outdriver = None
+    ds = None
     return
-
-
-# if __name__ == "__main__":
-#     input_image = r"E:\Test\sample_landsat5_image.tif"
-#     output_folder = r"E:\Test"
-#     lmax_list = [193.0, 365.0, 264.0, 221.0, 30.2, 16.5]
-#     lmin_list = [-1.520, -2.840, -1.170, -1.510, -0.370, -0.150]
-#     QCALMIN = 1
-#     QCALMAX = 255
-#     dn_to_radiance(input_image, output_folder, lmax_list, lmin_list, QCALMIN, QCALMAX)

src/raster_ops/components/export_land_class.py

@@ -36,10 +36,4 @@ def extract_lulc_class(lulc_raster, out_raster, class_number=10):
     result.SetGeoTransform(geot)
     result.SetProjection(geop)
     driver.CreateCopy(out_raster, result)
-    return "Process Completed!"
-
-
-# if __name__ == "__main__":
-#     lulc_raster = r"D:\MyWorkingFolder\sample_lulc.img"
-#     out_raster = r"D:\MyWorkingFolder\sample_forest.img"
-#     extract_lulc_class(lulc_raster, out_raster, class_number=10)
+    return

src/raster_ops/config.py

@@ -0,0 +1,28 @@
+import os
+
+from pathlib import Path
+from consts import *
+from raster_ops.config_entity import *
+
+
+class ConfigManager:
+    def __init__(self):
+        pass
+
+    def get_rasterdata_config(self):
+        rasterdata_config = RasterdataConfig(
+            raster_infile=RASTER_INPUT_FILE,
+            polygon_infile=CLIP_BOUNDARY_FILE,
+            clip_raster_file=CLIP_RASTER_OUTPUT_FILE,
+            out_point_file=OUT_POINT_FILE,
+            output_folder=RASTER_WORKING_DIR,
+            lmax_list=BANDS_LMAX_VALUES,
+            lmin_list=BANDS_LMIN_VALUES,
+            qcal_min=QCAL_MIN_VALUE,
+            qcal_max=QCAL_MAX_VALUE,
+            prefix=OUT_NAME_PREFIX,
+            lulc_raster_file=LULC_RASTER_FILE,
+            lulc_out_raster_file=LULC_OUT_CLASS_FILE,
+            lilc_class_code=LULC_CLASS_CODE
+        )
+        return rasterdata_config

src/raster_ops/config_entity.py

@@ -0,0 +1,19 @@
+from pathlib import Path
+from dataclasses import dataclass
+
+
+@dataclass(frozen=True)
+class RasterdataConfig:
+    raster_infile: Path
+    polygon_infile: Path
+    clip_raster_file: Path
+    out_point_file: Path
+    output_folder: Path
+    lmax_list: list
+    lmin_list: list
+    qcal_min: float
+    qcal_max: float
+    prefix: str
+    lulc_raster_file: Path
+    lulc_out_raster_file: Path
+    lilc_class_code: int

src/raster_ops/pipeline/clip_raster_pipeline.py

@@ -0,0 +1,43 @@
+from osgeo import gdal
+from logger import logger
+from raster_ops.config import ConfigManager
+from raster_ops.components.clip_raster import ClipRaster
+from utils import write_geotiff_file, np2gdal_dtype
+
+
+STAGE_NAME = "Clip raster by clip boundary extent"
+
+
+class ClipRasterPipeline:
+    def __init__(self) -> None:
+        pass
+
+    def main(self):
+        config = ConfigManager()
+        rstconfig = config.get_rasterdata_config()
+        obj = ClipRaster(rstconfig=rstconfig)
+        clip = obj.subset_by_extent()
+        # save raster
+        ds = gdal.Open(rstconfig.raster_infile, gdal.GA_ReadOnly)
+        src_geotransform = ds.GetGeoTransform()
+        src_projection = ds.GetProjectionRef()
+        ds = None
+        write_geotiff_file(
+            clip,
+            gt=src_geotransform,
+            sr=src_projection,
+            outfile_path=rstconfig.clip_raster_file,
+            dtype=np2gdal_dtype(str(clip.dtype))
+        )
+        return
+
+
+if __name__ == "__main__":
+    try:
+        logger.info(f">>>>> {STAGE_NAME} started <<<<<")
+        obj = ClipRasterPipeline()
+        obj.main()
+        logger.info(f">>>>> {STAGE_NAME} completed <<<<<\n")
+    except Exception as e:
+        logger.exception(e)
+        raise e

src/raster_ops/pipeline/dn2radiance_pipeline.py

@@ -0,0 +1,36 @@
+from logger import logger
+from raster_ops.config import ConfigManager
+from raster_ops.components.dn_to_radiance import dn_to_radiance
+
+
+STAGE_NAME = "Convert DN to Radiance of Satellite image bands"
+
+
+class DNToRadiancePipeline:
+    def __init__(self) -> None:
+        pass
+
+    def main(self):
+        config = ConfigManager()
+        rstconfig = config.get_rasterdata_config()
+        dn_to_radiance(
+            rstconfig.raster_infile,
+            rstconfig.output_folder,
+            rstconfig.lmax_list,
+            rstconfig.lmin_list,
+            rstconfig.qcal_min,
+            rstconfig.qcal_max,
+            prefix=rstconfig.prefix
+        )
+        return
+
+
+if __name__ == "__main__":
+    try:
+        logger.info(f">>>>> {STAGE_NAME} started <<<<<")
+        obj = DNToRadiancePipeline()
+        obj.main()
+        logger.info(f">>>>> {STAGE_NAME} completed <<<<<\n")
+    except Exception as e:
+        logger.exception(e)
+        raise e

src/raster_ops/pipeline/export_lulc_class_pipeline.py

@@ -0,0 +1,32 @@
+from logger import logger
+from raster_ops.config import ConfigManager
+from raster_ops.components.export_land_class import extract_lulc_class
+
+
+STAGE_NAME = "Export LULC individual class"
+
+
+class ExportLULCClassPipeline:
+    def __init__(self) -> None:
+        pass
+
+    def main(self):
+        config = ConfigManager()
+        rstconfig = config.get_rasterdata_config()
+        extract_lulc_class(
+            rstconfig.lulc_raster_file,
+            rstconfig.lulc_out_raster_file,
+            class_number=rstconfig.lilc_class_code
+        )
+        return
+
+
+if __name__ == "__main__":
+    try:
+        logger.info(f">>>>> {STAGE_NAME} started <<<<<")
+        obj = ExportLULCClassPipeline()
+        obj.main()
+        logger.info(f">>>>> {STAGE_NAME} completed <<<<<\n")
+    except Exception as e:
+        logger.exception(e)
+        raise e

src/raster_ops/pipeline/pixel2point_pipeline.py

@@ -0,0 +1,31 @@
+from logger import logger
+from raster_ops.components.convert_raster_to_point import raster_to_point
+from raster_ops.config import ConfigManager
+
+
+STAGE_NAME = "Convert gray-scale raster to point shapefile"
+
+
+class RasterToPointPipeline:
+    def __init__(self) -> None:
+        pass
+
+    def main(self):
+        config = ConfigManager()
+        rstconfig = config.get_rasterdata_config()
+        raster_to_point(
+            rstconfig.raster_infile,
+            rstconfig.out_point_file
+        )
+        return
+
+
+if __name__ == "__main__":
+    try:
+        logger.info(f">>>>> {STAGE_NAME} started <<<<<")
+        obj = RasterToPointPipeline()
+        obj.main()
+        logger.info(f">>>>> {STAGE_NAME} completed <<<<<\n")
+    except Exception as e:
+        logger.exception(e)
+        raise e

src/ref_scripts.py

@@ -340,3 +340,4 @@ def extract_pixel_value(imgfile, point_list):
 # http://scipy-lectures.github.io/advanced/image_processing/
 # http://stackoverflow.com/questions/4842613/merge-lists-that-share-common-elements
 # https://gis.stackexchange.com/questions/250077/overlapping-rasters-as-numpy-arrays
+# gdal raster datatype: http://ikcest-drr.osgeo.cn/tutorial/k8023