Installed rasterio through conda:

conda create -n gpkg-test -c conda-forge python=3.7 rasterio

This didn't work unfortunately Sean. No errors .... but nothing else
was wrote out as I was hoping. If you get a few minutes would you mind
testing to ensure I'm not going crazy? Thanks, James

On Thu, 27 Jun 2019 at 23:05, Sean Gillies via Groups.Io
<sean=mapbox.com@groups.io> wrote:
>
> Sorry, James, I messed up and used the GDAL vernacular. Replace YES with True and I think you'll be good to go.
>
> On Thu, Jun 27, 2019 at 3:09 PM James David Smith <james.david.smith@...> wrote:
>>
>> Thanks for the tip Sean, but that doesn't work for me.
>>
>> with rasterio.Env(TIFF_USE_OVR=True, GDAL_PAM_ENABLED=YES, ESRI_XML_PAM=YES):
>> ----- name 'YES' is not defined
>>
>> If I change 'YES' to True then the code runs .... but I'm not sure
>> it's recognizing the commands as the aux.xml file isn't created (but
>> the .ovr is).
>>
>> Anyone got a good idea?
>>
>>
>> On Thu, 27 Jun 2019 at 18:53, Sean Gillies via Groups.Io
>> <sean=mapbox.com@groups.io> wrote:
>> >
>> > Hi James,
>> >
>> > I have tried this myself, but I think you'll want to add both GDAL_PAM_ENABLED=YES and ESRI_XML_PAM=YES to the Env.
>> >
>> > https://trac.osgeo.org/gdal/wiki/ConfigOptions#GDAL_PAM_ENABLED
>> > https://gdal.org/drivers/raster/gtiff.html#configuration-options
>> >
>> >
>> > On Thu, Jun 27, 2019 at 11:11 AM James David Smith <james.david.smith@...> wrote:
>> >>
>> >> Hello,
>> >>
>> >> I'm trying to use rasterio to manipulate a file and then write it out,
>> >> along with accompanying ESRI-type files. So the output files will be:
>> >>
>> >> raster.tif,
>> >> raster.tif.aux.xml
>> >> raster.tif.ovr
>> >>
>> >> I think before writing the file I can do something like this to create
>> >> the ovr file.
>> >>
>> >> rasterio.Env(TIFF_USE_OVR=True)
>> >>
>> >> Is there a similar command to create the aux.xml file please?
>> >>
>> >> Thanks, James
>> >>
>> >>
>> >>
>> >
>> >
>> > --
>> > Sean Gillies
>> >
>>
>>
>>
>
>
> --
> Sean Gillies
>

I'm trying to add a raster layer to an existing GPKG with rasterio:

properties = {
        'driver':'GPKG',
        'width':array.shape[1],
        'height':array.shape[0],
        'count':1,
        'crs':crs,
        'transform':transform,
        'dtype':np.float32,
        'nodata':-9999,
        'layer':'rasterlayer'}
with rasterio.open('testing.gpkg', 'r+', **properties) as dst_dataset:
    dst_dataset.write(array)

However I get the following error:

Traceback (most recent call last):
  File "rasterio/_base.pyx", line 76, in rasterio._base.get_dataset_driver
  File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_OpenFailedError: 'testing.gpkg' not recognized as a supported file format.

I can do this with "gdal_translate" so I at least know the GPKG and rasters are all properly formatted:

gdal_translate -of GPKG raster.tiff testing.gpkg -co APPEND_SUBDATASET=YES -co RASTER_TABLE=rasterlayer


Side Question: Is there a way to format code in groups.io? I can't seem to find anything on the edit toolbar :-/

Thanks for the tip Sean, but that doesn't work for me.

with rasterio.Env(TIFF_USE_OVR=True, GDAL_PAM_ENABLED=YES, ESRI_XML_PAM=YES):
----- name 'YES' is not defined

If I change 'YES' to True then the code runs .... but I'm not sure
it's recognizing the commands as the aux.xml file isn't created (but
the .ovr is).

Anyone got a good idea?


On Thu, 27 Jun 2019 at 18:53, Sean Gillies via Groups.Io
<sean=mapbox.com@groups.io> wrote:
>
> Hi James,
>
> I have tried this myself, but I think you'll want to add both GDAL_PAM_ENABLED=YES and ESRI_XML_PAM=YES to the Env.
>
> https://trac.osgeo.org/gdal/wiki/ConfigOptions#GDAL_PAM_ENABLED
> https://gdal.org/drivers/raster/gtiff.html#configuration-options
>
>
> On Thu, Jun 27, 2019 at 11:11 AM James David Smith <james.david.smith@...> wrote:
>>
>> Hello,
>>
>> I'm trying to use rasterio to manipulate a file and then write it out,
>> along with accompanying ESRI-type files. So the output files will be:
>>
>> raster.tif,
>> raster.tif.aux.xml
>> raster.tif.ovr
>>
>> I think before writing the file I can do something like this to create
>> the ovr file.
>>
>> rasterio.Env(TIFF_USE_OVR=True)
>>
>> Is there a similar command to create the aux.xml file please?
>>
>> Thanks, James
>>
>>
>>
>
>
> --
> Sean Gillies
>

Hello,

I'm trying to use rasterio to manipulate a file and then write it out,
along with accompanying ESRI-type files. So the output files will be:

raster.tif,
raster.tif.aux.xml
raster.tif.ovr

I think before writing the file I can do something like this to create
the ovr file.

rasterio.Env(TIFF_USE_OVR=True)

Is there a similar command to create the aux.xml file please?

Thanks, James

OK, sorry, forget what I said - the same map tile that earlier was failing is now served fine. :/

As I mentioned, the canonical way to do this is to leverage functools.partial.

For example:

from functools import partial

from concurrent.futures import ProcessPoolExecutor

def worker_function(x, arg1, arg2):

    ...

if __name__ == '__main__':

    do_work = partial(worker_function, arg1=0, arg2='foo')

    with ProcessPoolExecutor(4) as executor:

        executor.map(do_work, [1, 2, 3])

The arguments have to be picklable, I think, but they should be in your case.

On 13/06/2019 16.29, javier lopatin via Groups.Io wrote:

Hi Dion, Thanks for the tip. I thought that that might be the problem. Any ideas of how to pass a function to executor.map that need several arguments? so far I'm only able to pass functions with a dstack as only argument, like:

def main(x):
    ...
    return something
executor.map(main, x)

so if I put the function main outside in the global scope I now need to have several args in the function:

def main(x, arg1, arg2, arg3, ...):
    ...
   return something

 

I have no idea how to pass these args to executor.map(). Any ideas?

I really appreciate this!

Cheers,

Javier

El jue., 13 jun. 2019 a las 15:12, Dion Häfner (<dion.haefner@...>) escribió:

Hey Javier,

this is not really a rasterio problem, mostly a Python peculiarity.

In general, you cannot pass functions defined in a local closure to a different process. Just move your main function to global scope and you should be fine. You can use functools.partial to fix some optional arguments of your main function and achieve the same thing as with your closure.

Best,
Dion

On 13/06/2019 15.05, javier lopatin via Groups.Io wrote:

Hi Guys!

I'm still working with the parallelization of raster processes with rasterio, and I came across a problem that you might know about. I've some raster stacks of time series data from which I want to obtain a interpolation signature. In practice I've a function to process the raster chunks in parallel with rasterio (calling a function main()):

def _parallel_process(inData, outData, main, count,  n_jobs=4, chuckSize=256):
    with rasterio.Env():
        with rasterio.open(inData) as src:
            profile = src.profile
            profile.update(blockxsize=chuckSize, blockysize=chuckSize,
                           count=count, dtype='float64', tiled=True)
            with rasterio.open(outData, "w", **profile) as dst:
                windows = [window for ij, window in dst.block_windows()]
                data_gen = (src.read(window=window) for window in windows)
                with concurrent.futures.ProcessPoolExecutor(
                    max_workers=n_jobs
                ) as executor:
                    for window, result in zip(
                        tqdm(windows), executor.map(main, data_gen)
                    ):
                        dst.write(result, window=window)

And I have another function that calls _parallel_process:

def PhenoShape(inData, outData, dates=None, nan_replace=None, rollWindow=None, nGS=46, chuckSize=256, n_jobs=4):
    def main(dstack):
        # load raster as a xarray
        xarray = xr.DataArray(dstack)
        xarray.coords['dim_0'] = dates
        xarray.coords['doy'] = xarray.dim_0.dt.dayofyear
        # sort basds according to day-of-the-year
        xarray = xarray.sortby('doy')
        # rearrange time dimension for smoothing and interpolation
        xarray['time'] = xarray['doy']
        # turn a value to NaN
        if nan_replace is not None:
            xarray = xarray.where(xarray.values != nan_replace)
        # rolling average using moving window
        if rollWindow is not None:
            xarray = xarray.rolling(dim_0=rollWindow, center=True).mean()
        # prepare inputs to getPheno
        x = xarray.doy.values
        y = xarray.values
        # get phenology shape accross the time axis
        return np.apply_along_axis(_getPheno, 0, y, x, nGS)
    # apply PhenoShape with parallel processing
    try:
        _parallel_process(inData, outData, main, nGS, n_jobs, chuckSize)
    except AttributeError:
        print('ERROR somehow.. :(')

 

This works really fine! But when I import the script as a library (functions inside phenology.py):

import phenology as phen
phen.PhenoShape(inData=inData, outData=outData, dates=dates, rollWindow=3,
    nan_replace=-32767, nGS=46, chuckSize=256, n_jobs=4)

I get the following error:

  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/queues.py", line 234, in _feed
    obj = _ForkingPickler.dumps(obj)
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/reduction.py", line 51, in dumps
    cls(buf, protocol).dump(obj)
AttributeError: Can't pickle local object 'PhenoShape.<locals>.main'
Traceback (most recent call last):
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/queues.py", line 234, in _feed
    obj = _ForkingPickler.dumps(obj)
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/reduction.py", line 51, in dumps
    cls(buf, protocol).dump(obj)
AttributeError: Can't pickle local object 'PhenoShape.<locals>.main'
Traceback (most recent call last):
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/queues.py", line 234, in _feed
    obj = _ForkingPickler.dumps(obj)
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/reduction.py", line 51, in dumps
    cls(buf, protocol).dump(obj)
AttributeError: Can't pickle local object 'PhenoShape.<locals>.main'
Traceback (most recent call last):
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/queues.py", line 234, in _feed
    obj = _ForkingPickler.dumps(obj)
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/reduction.py", line 51, in dumps
    cls(buf, protocol).dump(obj)
AttributeError: Can't pickle local object 'PhenoShape.<locals>.main'
Traceback (most recent call last):
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/queues.py", line 234, in _feed
    obj = _ForkingPickler.dumps(obj)
  File "/home/javier/miniconda3/lib/python3.6/multiprocessing/reduction.py", line 51, in dumps
    cls(buf, protocol).dump(obj)
AttributeError: Can't pickle local object 'PhenoShape.<locals>.main'

Any ideas? The main idea here is to call this functions to have cleaner script at the end.

Many thanks in advance.

Cheers,

Javier L.

El jue., 23 may. 2019 a las 11:13, javier lopatin via Groups.Io (<javierlopatin=gmail.com@groups.io>) escribió:

Hi Sean, thanks a lot. This was very helpfull indeed, and the performance improved a bit when changing the size to 512. 

Grate job on the library, congrats!

Cheers,

Javier

El mié., 22 may. 2019 a las 23:36, Alan Snow (<alansnow21@...>) escribió:

You are indeed correct Sean. Thanks for catching that and providing the correct answer! It seems I skimmed through the original message too fast.

--

Javier Lopatin
Institute of Geography and Geoecology
Karlsruhe Institute of Technology (KIT)
javier.lopatin@...

--

Javier Lopatin
Institute of Geography and Geoecology
Karlsruhe Institute of Technology (KIT)
javier.lopatin@...

--

Javier Lopatin
Institute of Geography and Geoecology
Karlsruhe Institute of Technology (KIT)
javier.lopatin@...