Revised: Don't require the seek call:

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

For reference, doing:

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    profile['driver'] = 'GTiff'
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

results in a correctly geocoded GeoTiff. Appears it is an issue with .jp2s and Profiles/MemoryFiles.

On further inspection, when I switched to saving my outputs as GeoTIFF, the resultant images are both geo-coded fine (Disk vs. In Memory to S3).

I'm guessing something is going wrong when .jp2 is involved...

Any suggestions on next steps?

.jp2 is both preferable (because of file sizes) and a requirement for the work I'm currently doing...

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    stack.seek(0)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

If I download that image and open it in QGIS, the image has no geo-coding at all.

If I process the same image, but write to disk like:

profile = source_image.profile.copy()

with rasterio.open('thing.jp2', 'w', **profile) as stack_out:
    stack_out.write(source_image.read(1), 1))

I get a geo-coded image and it appears where I'd expect.

If I compare file sizes:
MemoryFile(): 49,211KB
Disk: 49,214KB

I can only assume I'm losing any profile information when using a MemoryFile, can someone direct me to what I'm doing wrong?

Thanks!

On dimanche 26 avril 2020 03:04:56 CEST carlogarro@... wrote:

> Hello, I need to read many huge datasets and the speed time is very

> important to avoid a bottleneck. I have to read a tiff file that has 20

> bands, and a window of 224,224. Now I am doing like this, and it takes

> approx 0.8seconds.

>

> with rasterio.open('./sentinel.tif') as src:

> sentinel1_1 = src.read(window=window)

>

> What I realized is that if I try to read only one of the bands the required

> time is approx the same, but when reading a tiff of only one band the

> amount of time is 10 times shorter.

>

> Can I do something to speed it up? Maybe read bands in parallel, I don't

> really know.

If you have control on how the creation of the TIFF file, make sure it uses Band interleaving instead of Pixel interleaving

For example, with gdal_translate can be done with -co INTERLEAVE=BAND

If the file is not tiled, adding tiling might also help.

I see in https://rasterio.readthedocs.io/en/latest/topics/profiles.html a pure rasterio way of creating such file

--

Spatialys - Geospatial professional services

http://www.spatialys.com

You can try the following:

So my issue is resolved. The issue was caused by having a GDAL installation on my windows machine, as well as all the GDAL paths set up. This was conflicting with the GDAL setup in each conda environment.

I uninstalled all GDAL core and bindings from Windows itself, and then only installed GDAL in the Anaconda environment. No system variables were set after the GDAL was installed in the environment.

On my windows PC, I am running Anaconda to handle my environments.

My guess is that the issue resides in the GDAL installation on windows, that was performed prior to installing anaconda.
PROJ_LIB was set in the system variable and pointed to the proj.db file in the GDAL folder.
If I create a new environment in Conda, and then install rasterio / gdal / geopandas, it is likely to pick up the GDAL installation currently on windows, and assume that PROJ_LIB is correct (and proj.db)? Is this likely to be leading to the error?

In terms of Conda environments and GDAL/Rasterio, should we treat each environment with its own GDAL installation and own routes to GDAL_DATA path's?

I appreciate the assistance, but how can I ensure correct environment mappings of PROJ to the proj.db that up to date with the correct db schema in the proj.db file? Is there a good rule of thumb to not stuff this up in future?

> ERROR 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326:

> SQLite error on SELECT name, ellipsoid_auth_name, ellipsoid_code,

> prime_meridian_auth_name, prime_meridian_code, area_of_use_auth_name,

> area_of_use_code, publication_date, deprecated FROM geodetic_datum WHERE

> auth_name = ? AND code = ?: no such column: publication_date Traceback

Smells like you have PROJ 6.x.y version pointing to the proj.db of another PROJ 6.z.t version, since there were a few database schema changes. Make sure your PROJ_LIB is correctly set.

--

Spatialys - Geospatial professional services

http://www.spatialys.com

       0.0, -15.0, 1438430.0), 'tiled': False, 'interleave': 'pixel'}

Hi everyone. I started using rasterio a few days ago to read some DSM/DTM and process them using numba.

I was trying to debug some weird behaviour in my numba function and i realized that the "read" function in Rasterio actually flips the axis of my dataset.
The file I am opening is a Geotiff, with EPSG:3003 projection.

This how i reproduced this behaviour:

dataset = rasterio.open("raster.tif")
band = dataset.read(1)
print(dataset.width)    -> 84
print(dataset.heigth)  -> 44
print(band.shape)   -> (44,84)


Is this behaviour normal? I don't understand why my data are being flipped in this way.

Thanks, Gabriel