Error reprojecting using gdal

[Katospiegel]

Hello,

Thanks for making your work open source. I saw the results, and they are amazing. I'm trying to do a little test on Mac using gdal 3.5.1. However, I have the error below. Could it be an error related to the gdal version?

Thanks a lot in advance,
Carlos

2022-07-26 23:44:22,796 - SIAC-V2.3.6 - INFO - Preprocessing for S2B_MSIL1C_20181225T143749_N0207_R096_T19HCC_20181225T175914.SAFE
2022-07-26 23:44:22,797 - SIAC-V2.3.6 - INFO - Doing per pixel angle resampling.
2022-07-26 23:44:58,394 - SIAC-V2.3.6 - INFO - Getting cloud mask.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
File <timed exec>:4, in <module>

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/SIAC/SIAC_S2.py:44, in SIAC_S2(s2_t, send_back, mcd43, vrt_dir, aoi, global_dem, cams_dir, jasmin, Gee)
     27 def SIAC_S2(s2_t, send_back = False, mcd43 = home + '/MCD43/', vrt_dir = home + '/MCD43_VRT/', aoi = None, 
     28              global_dem  = None, cams_dir = None, jasmin = False, Gee = True):
     29     '''
     30     if not os.path.exists(file_path + '/emus/'):
     31         os.mkdir(file_path + '/emus/')
   (...)
     42         parmap(f, to_down)
     43     '''
---> 44     rets = s2_pre_processing(s2_t, cams_dir, global_dem)
     45     aero_atmos = []
     46     for ret in rets:

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/SIAC/s2_preprocessing.py:226, in s2_pre_processing(s2_dir, cams_dir, dem)
    223 off = off[[0,1,3,4,7,8,9,10,11,12], None, None]# (band, nx, ny)
    225 logger.info('Getting cloud mask.')
--> 226 cloud = do_cloud(cloud_bands, cloud_name, ref_scale=scale, ref_offset=off)
    227 cloud_mask = cloud > 0.6# 0.4 is the sentinel hub default
    228 clean_pixel = ((cloud>=0).sum() - cloud_mask.sum()) / 3348900. * 100.

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/SIAC/s2_preprocessing.py:44, in do_cloud(cloud_bands, cloud_name, ref_scale, ref_offset)
     40 def do_cloud(cloud_bands, cloud_name = None, ref_scale = 1/10000., ref_offset = 0.):
     42     cl = Booster(model_file=model_filename)
---> 44     toas = [reproject_data(str(band), cloud_bands[0], dstNodata=0, resample=5).data for band in cloud_bands]
     45     # print(ref_scale, ref_offset)
     46     toas = np.array(toas) * ref_scale + ref_offset

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/SIAC/s2_preprocessing.py:44, in <listcomp>(.0)
     40 def do_cloud(cloud_bands, cloud_name = None, ref_scale = 1/10000., ref_offset = 0.):
     42     cl = Booster(model_file=model_filename)
---> 44     toas = [reproject_data(str(band), cloud_bands[0], dstNodata=0, resample=5).data for band in cloud_bands]
     45     # print(ref_scale, ref_offset)
     46     toas = np.array(toas) * ref_scale + ref_offset

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/SIAC/reproject.py:79, in reproject_data.__init__(self, source_img, target_img, dstSRS, srcNodata, dstNodata, outputType, verbose, xmin, xmax, ymin, ymax, xRes, yRes, xSize, ySize, resample)
     77     raster_wkt = g.GetProjection()
     78     dstSRS.ImportFromWkt(raster_wkt)
---> 79     self.g = gdal.Warp('', self.source_img, format = 'MEM', outputBounds = [xmin, ymin, xmax, ymax], dstNodata=self.dstNodata, warpOptions = ['NUM_THREADS=ALL_CPUS'],\
     80                         xRes = self.xRes, yRes = self.yRes, dstSRS = dstSRS, outputType = self.outputType, srcNodata = self.srcNodata, resampleAlg = self.resample)
     82 else:
     83     self.g = gdal.Warp('', self.source_img, format = 'MEM', outputBounds = [self.xmin, self.ymin, \
     84                        self.xmax, self.ymax], xRes = self.xRes, yRes = self.yRes, dstSRS = self.dstSRS, warpOptions = ['NUM_THREADS=ALL_CPUS'],\
     85                        copyMetadata=True, outputType = self.outputType, dstNodata=self.dstNodata, srcNodata = self.srcNodata, resampleAlg = self.resample)

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/osgeo/gdal.py:660, in Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)
    650 """ Warp one or several datasets.
    651     Arguments are :
    652       destNameOrDestDS --- Output dataset name or object
   (...)
    656       other keywords arguments of gdal.WarpOptions()
    657     If options is provided as a gdal.WarpOptions() object, other keywords are ignored. """
    659 if 'options' not in kwargs or isinstance(kwargs['options'], (list, str)):
--> 660     (opts, callback, callback_data) = WarpOptions(**kwargs)
    661 else:
    662     (opts, callback, callback_data) = kwargs['options']

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/osgeo/gdal.py:647, in WarpOptions(options, format, outputBounds, outputBoundsSRS, xRes, yRes, targetAlignedPixels, width, height, srcSRS, dstSRS, coordinateOperation, srcAlpha, dstAlpha, warpOptions, errorThreshold, warpMemoryLimit, creationOptions, outputType, workingType, resampleAlg, srcNodata, dstNodata, multithread, tps, rpc, geoloc, polynomialOrder, transformerOptions, cutlineDSName, cutlineLayer, cutlineWhere, cutlineSQL, cutlineBlend, cropToCutline, copyMetadata, metadataConflictValue, setColorInterpretation, overviewLevel, callback, callback_data)
    644 if return_option_list:
    645     return new_options
--> 647 return (GDALWarpAppOptions(new_options), callback, callback_data)

File ~/anaconda3/envs/geospatial_3.8/lib/python3.8/site-packages/osgeo/gdal.py:4318, in GDALWarpAppOptions.__init__(self, *args)
   4316 def __init__(self, *args):
   4317     r"""__init__(GDALWarpAppOptions self, char ** options) -> GDALWarpAppOptions"""
-> 4318     _gdal.GDALWarpAppOptions_swiginit(self, _gdal.new_GDALWarpAppOptions(*args))

RuntimeError: Translating source or target SRS failed:

[MarcYin]

Hi,

Thanks a lot for your interests on SIAC.

I have tested with your file and it seems working for me:

from SIAC import SIAC_S2
global_dem = '/vsicurl/https://gws-access.jasmin.ac.uk/public/nceo_ard/DEM_V3/global_dem.vrt'
cams_dir = '/vsicurl/https://gws-access.jasmin.ac.uk/public/nceo_ard/cams/'
SIAC_S2('./Downloads/S2B_MSIL1C_20181225T143749_N0207_R096_T19HCC_20181225T175914.SAFE/', global_dem = global_dem, cams_dir = cams_dir)
2022-07-27 12:35:08,713 - SIAC-V2.3.6 - INFO - Preprocessing for S2B_MSIL1C_20181225T143749_N0207_R096_T19HCC_20181225T175914.SAFE
2022-07-27 12:35:08,713 - SIAC-V2.3.6 - INFO - Doing per pixel angle resampling.

2022-07-27 12:35:36,464 - SIAC-V2.3.6 - INFO - Getting cloud mask.
2022-07-27 12:36:11,832 - SIAC-V2.3.6 - INFO - Valid pixel percentage: 100.00
2022-07-27 12:36:11,832 - SIAC-V2.3.6 - INFO - Clean pixel percentage: 99.94
2022-07-27 12:36:11,993 - SIAC-V2.3.6 - INFO - Starting atmospheric corretion for S2B_MSIL1C_20181225T143749_N0207_R096_T19HCC_20181225T175914.SAFE
2022-07-27 12:36:11,997 - SIAC-V2.3.6 - INFO - Getting MCD43 from GEE

I have gdal and libgdal version 3.5.0, which should be pretty much the same as yours. Could you please firstly check the gdal package is working properly?

Try gdalwarp with:

gdalwarp /vsicurl/https://gws-access.jasmin.ac.uk/public/nceo_ard/DEM_V3/global_dem.vrt out_dem.tif -te 300000 6190240 409800 6300040 -t_srs EPSG:32719

Thanks,

Marc.

[Katospiegel]

Hello Marc,

Thanks a lot for the quick answer. Yes, the gdalwarp is working with that command.

(base) kato@nostromo:~/pro/wekeo/git/portrait_of_a_lakes_death$ gdalwarp /vsicurl/https://gws-access.jasmin.ac.uk/public/nceo_ard/DEM_V3/global_dem.vrt out_dem.tif -te 300000 6190240 409800 6300040 -t_srs EPSG:32719
Creating output file that is 3063P x 3063L.
Processing /vsicurl/https://gws-access.jasmin.ac.uk/public/nceo_ard/DEM_V3/global_dem.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

Could it be that the error is arising upstream? I'm going to try to print some logs to what is given exactly what python is giving to the wrapper.

Thanks,
Carlos