search

centreborelli / tsd

Time Series Downloader for satellite images
GNU Affero General Public License v3.0
40 stars 19 forks source link

Replace gdal with rasterio #20

Closed glostis closed 4 years ago

glostis commented 4 years ago

The main goal of this PR is to remove tsd's dependency on GDAL.

The only remaining function using GDAL was tsd.utils.crop_with_gdal_translate. It has been replaced by tsd.utils.rasterio_geo_crop which uses rasterio to obtain the same result.

Another refactoring that was done was to call the new cropping function with an epsg argument instead of utm_zone and lat_band.

Also, the function doesn't need to specifically handle gs:// input paths, thanks to #17.

You can test that the new function gives (nearly) the same result as the previous one, by running:

  1. On this branch:
Expand test code ```python from tsd.utils import rasterio_geo_crop rasterio_geo_crop( outpath="tests/l8_s3_rio.tif", inpath="https://s3-us-west-2.amazonaws.com/landsat-pds/c1/L8/176/039/LC08_L1TP_176039_20180411_20180417_01_T1/LC08_L1TP_176039_20180411_20180417_01_T1_B8.TIF", ulx=318753.0277473631, uly=3319235.8866136344, lrx=321313.0361352647, lry=3316675.838327111, epsg=32636, output_type=None, ) rasterio_geo_crop( outpath="tests/l8_gcp_rio.tif", inpath="https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/176/039/LC08_L1TP_176039_20180411_20180417_01_T1/LC08_L1TP_176039_20180411_20180417_01_T1_B8.TIF", ulx=318753.0277473631, uly=3319235.8866136344, lrx=321313.0361352647, lry=3316675.838327111, epsg=32636, output_type=None, ) rasterio_geo_crop( outpath="tests/s2_s3_rio.tif", inpath="s3://sentinel-s2-l1c/tiles/36/R/UU/2018/4/9/0/B04.jp2", ulx=318780.0, uly=3319260.0, lrx=321300.0, lry=3316680.0, epsg=32636, output_type="UInt16", ) rasterio_geo_crop( outpath="tests/s2_gcp_rio.tif", inpath="https://storage.googleapis.com/gcp-public-data-sentinel-2/tiles/36/R/UU/S2A_MSIL1C_20180409T082601_N0206_R021_T36RUU_20180409T121946.SAFE/GRANULE/L1C_T36RUU_A014604_20180409T082924/IMG_DATA/T36RUU_20180409T082601_B04.jp2", ulx=318780.0, uly=3319260.0, lrx=321300.0, lry=3316680.0, epsg=32636, output_type="UInt16", ) ```
  1. On master:
    Expand test code
from tsd.utils import crop_with_gdal_translate

crop_with_gdal_translate(
    outpath="tests/l8_s3_gdal.tif",
    inpath="https://s3-us-west-2.amazonaws.com/landsat-pds/c1/L8/176/039/LC08_L1TP_176039_20180411_20180417_01_T1/LC08_L1TP_176039_20180411_20180417_01_T1_B8.TIF",
    ulx=318753.0277473631,
    uly=3319235.8866136344,
    lrx=321313.0361352647,
    lry=3316675.838327111,
    utm_zone=36,
    lat_band="R",
    output_type=None,
)

crop_with_gdal_translate(
    outpath="tests/l8_gcp_gdal.tif",
    inpath="https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/176/039/LC08_L1TP_176039_20180411_20180417_01_T1/LC08_L1TP_176039_20180411_20180417_01_T1_B8.TIF",
    ulx=318753.0277473631,
    uly=3319235.8866136344,
    lrx=321313.0361352647,
    lry=3316675.838327111,
    utm_zone=36,
    lat_band="R",
    output_type=None,
)

crop_with_gdal_translate(
    outpath="tests/s2_s3_gdal.tif",
    inpath="s3://sentinel-s2-l1c/tiles/36/R/UU/2018/4/9/0/B04.jp2",
    ulx=318780.0,
    uly=3319260.0,
    lrx=321300.0,
    lry=3316680.0,
    utm_zone=36,
    lat_band="R",
    output_type="UInt16",
)

crop_with_gdal_translate(
    outpath="tests/s2_gcp_gdal.tif",
    inpath="https://storage.googleapis.com/gcp-public-data-sentinel-2/tiles/36/R/UU/S2A_MSIL1C_20180409T082601_N0206_R021_T36RUU_20180409T121946.SAFE/GRANULE/L1C_T36RUU_A014604_20180409T082924/IMG_DATA/T36RUU_20180409T082601_B04.jp2",
    ulx=318780.0,
    uly=3319260.0,
    lrx=321300.0,
    lry=3316680.0,
    utm_zone=36,
    lat_band="R",
    output_type="UInt16",
)

Expand diff between results ``` diff <(gdalinfo tests/l8_gcp_gdal.tif) <(gdalinfo tests/l8_gcp_rio.tif) ``` ```diff 2c2 < Files: tests/l8_gcp_gdal.tif --- > Files: tests/l8_gcp_rio.tif 30c30 < AREA_OR_POINT=Point --- > AREA_OR_POINT=Area 40c40 < Band 1 Block=171x23 Type=UInt16, ColorInterp=Gray --- > Band 1 Block=256x256 Type=UInt16, ColorInterp=Gray ``` ``` diff <(gdalinfo tests/s2_gcp_gdal.tif) <(gdalinfo tests/s2_gcp_rio.tif) ``` ```diff 2c2 < Files: tests/s2_gcp_gdal.tif --- > Files: tests/s2_gcp_rio.tif 40c40 < Band 1 Block=252x16 Type=UInt16, ColorInterp=Gray --- > Band 1 Block=1024x1024 Type=UInt16, ColorInterp=Gray ```

Sorry about the readability of this PR... Maybe you'll find it easier to review the diff commit by commit.

glostis commented 4 years ago

Could be addressed in a future PR: clean up/remove files that refer to the installation process with GDAL, which is not needed anymore: