Date   

Re: Rasterio read outshape resample does not give correct array

Hobart, Geordie (NRCan/RNCan)
 

Hello.

I think that your problem is a little more complex than just resampling. The bilinear interpolation averages the surrounding cells to create a smooth image.

You need to sum the four surrounding cells, and populate the new grid with the result.

Not sure exactly how to do it but a clever google should find what you need.

Good Luck.

 

 

From: main@rasterio.groups.io <main@rasterio.groups.io> On Behalf Of tinak.contact@...
Sent: March 16, 2021 14:26
To: main@rasterio.groups.io
Subject: [rasterio] Rasterio read outshape resample does not give correct array

 

[Edited Message Follows]

Hi all, 

I have a population grid and want to decrease the resolution by a factor 2. So I would like to get a raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same overall number of population.

 

This should be a straightforward operation, but the population numbers change incorrectly if I follow the documentation. When the resampled array is written to file, the numbers are incorrect. I have a population grid and want to decrease the resolution by a factor 2. So I would like to get a raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same number of population.

 

Shouldn't the resample method yield an array with resampled cell values to get the same number of total population? Any tips and hints would be appreciated!

 

factor = 2
 
with rasterio.open("mwi_ppp_2020.tif") as dataset:
 
    arr = dataset.read(1,masked=True)
    #Show total population of Malawi 
    print(arr.sum())
 
    data = dataset.read(1, masked=True, resampling=Resampling.bilinear,
        out_shape=(
            int(dataset.height * factor),
            int(dataset.width * factor)))
 
     #Show total population of array 
    print(data.sum())
 
     # scale image transform
    transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[1]),

 

            (dataset.height / data.shape[0])
                                )


Rasterio read outshape resample does not give correct array

tinak.contact@...
 
Edited

Hi all, 

I have a population grid and want to decrease the resolution by a factor 2. So I would like to get a raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same overall number of population.
 
This should be a straightforward operation, but the population numbers change incorrectly if I follow the documentation. When the resampled array is written to file, the numbers are incorrect.  
 
Shouldn't the resample method yield an array with resampled cell values to get the same number of total population? Any tips and hints would be appreciated!

factor = 2

with rasterio.open("mwi_ppp_2020.tif") as dataset:

    arr = dataset.read(1,masked=True)
    #Show total population of Malawi 
    print(arr.sum())

    data = dataset.read(1, masked=True, resampling=Resampling.bilinear,
        out_shape=(
            int(dataset.height * factor),
            int(dataset.width * factor)))

     #Show total population of array 
    print(data.sum())

     # scale image transform
    transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[1]),


            (dataset.height / data.shape[0])
                                )


Re: cannot correctly change a TIFF's interleave

guohoucai@...
 

Thanks! It really helps me a lot. Hope you have a nice day :)


Re: cannot correctly change a TIFF's interleave

Sean Gillies
 

The "profile" attribute of a dataset is a dictionary that describes key characteristics of the dataset. You can use this dictionary in any way you like. Modifying it does not affect the dataset that it originated from. Common rasterio usage is to get the profile of one dataset and then use that as a baseline for creation of a new dataset. Here's an example: https://github.com/mapbox/rasterio/blob/master/rasterio/rio/convert.py#L62-L77.

Does this make sense? I hope it helps.


On Fri, Mar 5, 2021 at 6:46 PM <guohoucai@...> wrote:
Thanks for your help! It really helps me a lot! But I got another question: rasterio.shutil.copy() could only copy a raster to a new destination, what should I do if I wanna use a raster's profile for the creation of a new raster? Like I use a multi-bands TIFF file to calculate the NDVI and then save it, the multi-bands TIFF's profile is useful for the NDVI TIFF file cause they have the same "crs" and "transform". I'm a rasterio rookie and I appreciate your generosity a lot!


--
Sean Gillies


Re: cannot correctly change a TIFF's interleave

guohoucai@...
 

Thanks for your help! It really helps me a lot! But I got another question: rasterio.shutil.copy() could only copy a raster to a new destination, what should I do if I wanna use a raster's profile for the creation of a new raster? Like I use a multi-bands TIFF file to calculate the NDVI and then save it, the multi-bands TIFF's profile is useful for the NDVI TIFF file cause they have the same "crs" and "transform". I'm a rasterio rookie and I appreciate your generosity a lot!


Re: cannot correctly change a TIFF's interleave

Sean Gillies
 

Hello,

On Thu, Mar 4, 2021 at 7:24 PM <guohoucai@...> wrote:
I'm looking for some help about how to change a TIFF's interleave. My code could change a TIFF's nodata value but cannot succeed on the interleave, could anyone give me some advice? I'll appreciate your kindness.
_._,_._,_

The nodata value is nothing but an annotation of the data and is simple to change. To change the interleaving of a GeoTIFF, it is required to rewrite the file entirely. https://rasterio.readthedocs.io/en/latest/api/rasterio.shutil.html?highlight=copy#rasterio.shutil.copy is the function to use for this purpose. I'm sorry that the documentation of the shutil module is incomplete, we're working on it.

I recommend this usage:

with rasterio.open(test_img) as src:
    prof = src.profile
    prof.update(interleave="pixel", nodata=0)
    rasterio.shutil.copy(src, save_path, driver="GTiff", **prof)

--
Sean Gillies


cannot correctly change a TIFF's interleave

guohoucai@...
 

I'm looking for some help about how to change a TIFF's interleave. My code could change a TIFF's nodata value but cannot succeed on the interleave, could anyone give me some advice? I'll appreciate your kindness.


Rasterio 1.2.1

Sean Gillies
 

Hi all,

Rasterio 1.2.1 is on PyPI now. Big thanks to Ryan Grout, Denis Rykov, Kirill Kouzoubov, Martijn Visser, and Alan Snow.


--
Sean Gillies


Re: Clarification on usage with QGIS

Ari Meyer
 

Much appreciated, Christina.  Will try this out.

Cheers!
Ari


Re: rasterio merge producing strange results

Paolo Corti
 

Dear Sean
I had the chance to run one process with the above version, and it
looks like the bug went away :)
We are now looking forward to use the 1.2.1 release
Thanks a lot you and all the team for all the wonderful work
cheers
Paolo

On Tue, Mar 2, 2021 at 11:31 AM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:

The maint-1.2 branch, specifically commit b962f10, is the thing to try. If you're willing, it would be wonderful to hear how it goes.

On Tue, Mar 2, 2021 at 9:03 AM Paolo Corti <pcorti@gmail.com> wrote:

Thanks a lot Sean, it is great to know that.
If you like I could have a test run with my data. Is the code change
already committed? If so, which branch/commit should I use for the
purpose?
Have a great day, cheers
Paolo

On Tue, Mar 2, 2021 at 10:50 AM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:

Hi Paolo,

The fix (we believe) will be in version 1.2.1 https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L9. Later this week.

On Tue, Mar 2, 2021 at 8:45 AM Paolo Corti <pcorti@gmail.com> wrote:

Hi Denis
thanks a lot for helping. I have tried to use rasterio 1.2.0 but
unfortunately I am experiencing the same problem
Do you know if there is anything which can be done to avoid this
issue? Otherwise I will need to change my code and maybe use GDAL
virtual rasters to do the merge process
thanks in advance!
Paolo

On Tue, Mar 2, 2021 at 8:52 AM Denis Rykov <rykovd@gmail.com> wrote:

Hi Paolo. Can you reproduce the issue with the latest version of rasterio?

On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@gmail.com> wrote:

Hello all

I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).

I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.

Any idea why this is happening? Do you have any suggestions for fixing this?

thanks a lot in advance
Paolo



[1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1





--
Sean Gillies


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1





--
Sean Gillies


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Re: rasterio merge producing strange results

Sean Gillies
 

The maint-1.2 branch, specifically commit b962f10, is the thing to try. If you're willing, it would be wonderful to hear how it goes.

On Tue, Mar 2, 2021 at 9:03 AM Paolo Corti <pcorti@...> wrote:
Thanks a lot Sean, it is great to know that.
If you like I could have a test run with my data. Is the code change
already committed? If so, which branch/commit should I use for the
purpose?
Have a great day, cheers
Paolo

On Tue, Mar 2, 2021 at 10:50 AM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:
>
> Hi Paolo,
>
> The fix (we believe) will be in version 1.2.1 https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L9. Later this week.
>
> On Tue, Mar 2, 2021 at 8:45 AM Paolo Corti <pcorti@...> wrote:
>>
>> Hi Denis
>> thanks a lot for helping. I have tried to use rasterio 1.2.0 but
>> unfortunately I am experiencing the same problem
>> Do you know if there is anything which can be done to avoid this
>> issue? Otherwise I will need to change my code and maybe use GDAL
>> virtual rasters to do the merge process
>> thanks in advance!
>> Paolo
>>
>> On Tue, Mar 2, 2021 at 8:52 AM Denis Rykov <rykovd@...> wrote:
>> >
>> > Hi Paolo. Can you reproduce the issue with the latest version of rasterio?
>> >
>> > On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@...> wrote:
>> >>
>> >> Hello all
>> >>
>> >> I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).
>> >>
>> >> I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
>> >> This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.
>> >>
>> >> Any idea why this is happening? Do you have any suggestions for fixing this?
>> >>
>> >> thanks a lot in advance
>> >> Paolo
>> >>
>> >>
>> >>
>> >> [1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result
>> >>
>> >> --
>> >> Paolo Corti
>> >> Geospatial software developer
>> >> web: http://www.paolocorti.net
>> >> twitter: @capooti
>> >> skype: capooti
>> >> #drt3jc1
>> >
>> >
>>
>>
>>
>> --
>> Paolo Corti
>> Geospatial software developer
>> web: http://www.paolocorti.net
>> twitter: @capooti
>> skype: capooti
>> #drt3jc1
>>
>>
>>
>>
>>
>
>
> --
> Sean Gillies
>



--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1







--
Sean Gillies


Re: rasterio merge producing strange results

Paolo Corti
 

Thanks a lot Sean, it is great to know that.
If you like I could have a test run with my data. Is the code change
already committed? If so, which branch/commit should I use for the
purpose?
Have a great day, cheers
Paolo

On Tue, Mar 2, 2021 at 10:50 AM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:

Hi Paolo,

The fix (we believe) will be in version 1.2.1 https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L9. Later this week.

On Tue, Mar 2, 2021 at 8:45 AM Paolo Corti <pcorti@gmail.com> wrote:

Hi Denis
thanks a lot for helping. I have tried to use rasterio 1.2.0 but
unfortunately I am experiencing the same problem
Do you know if there is anything which can be done to avoid this
issue? Otherwise I will need to change my code and maybe use GDAL
virtual rasters to do the merge process
thanks in advance!
Paolo

On Tue, Mar 2, 2021 at 8:52 AM Denis Rykov <rykovd@gmail.com> wrote:

Hi Paolo. Can you reproduce the issue with the latest version of rasterio?

On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@gmail.com> wrote:

Hello all

I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).

I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.

Any idea why this is happening? Do you have any suggestions for fixing this?

thanks a lot in advance
Paolo



[1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1





--
Sean Gillies


--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Re: rasterio merge producing strange results

Sean Gillies
 

Hi Paolo,

The fix (we believe) will be in version 1.2.1 https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L9. Later this week.

On Tue, Mar 2, 2021 at 8:45 AM Paolo Corti <pcorti@...> wrote:
Hi Denis
thanks a lot for helping. I have tried to use rasterio 1.2.0 but
unfortunately I am experiencing the same problem
Do you know if there is anything which can be done to avoid this
issue? Otherwise I will need to change my code and maybe use GDAL
virtual rasters to do the merge process
thanks in advance!
Paolo

On Tue, Mar 2, 2021 at 8:52 AM Denis Rykov <rykovd@...> wrote:
>
> Hi Paolo. Can you reproduce the issue with the latest version of rasterio?
>
> On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@...> wrote:
>>
>> Hello all
>>
>> I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).
>>
>> I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
>> This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.
>>
>> Any idea why this is happening? Do you have any suggestions for fixing this?
>>
>> thanks a lot in advance
>> Paolo
>>
>>
>>
>> [1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result
>>
>> --
>> Paolo Corti
>> Geospatial software developer
>> web: http://www.paolocorti.net
>> twitter: @capooti
>> skype: capooti
>> #drt3jc1
>
>



--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1







--
Sean Gillies


Re: rasterio merge producing strange results

Paolo Corti
 

Hi Denis
thanks a lot for helping. I have tried to use rasterio 1.2.0 but
unfortunately I am experiencing the same problem
Do you know if there is anything which can be done to avoid this
issue? Otherwise I will need to change my code and maybe use GDAL
virtual rasters to do the merge process
thanks in advance!
Paolo

On Tue, Mar 2, 2021 at 8:52 AM Denis Rykov <rykovd@gmail.com> wrote:

Hi Paolo. Can you reproduce the issue with the latest version of rasterio?

On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@gmail.com> wrote:

Hello all

I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).

I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.

Any idea why this is happening? Do you have any suggestions for fixing this?

thanks a lot in advance
Paolo



[1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1
--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Re: rasterio merge producing strange results

Denis Rykov
 

Hi Paolo. Can you reproduce the issue with the latest version of rasterio? 


On Tue, Mar 2, 2021, 2:46 PM Paolo Corti <pcorti@...> wrote:
Hello all

I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).

I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.

Any idea why this is happening? Do you have any suggestions for fixing this?

thanks a lot in advance
Paolo

Screen Shot 2021-03-02 at 8.37.16 AM.png

[1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


rasterio merge producing strange results

Paolo Corti
 

Hello all

I am experiencing strange results with rasterio.merge module (rasterio==1.1.5).

I am merging a very large number of small images (in chunks of 100), and the output always has strange stripes in the central area of the output image (see attached file).
This is happening for several areas of the world. In some cases I have found that setting the precision parameter to 50, as suggested here [1], seems to fix the problem, however there are cases where I still get the stripes.

Any idea why this is happening? Do you have any suggestions for fixing this?

thanks a lot in advance
Paolo

Screen Shot 2021-03-02 at 8.37.16 AM.png

[1] https://gis.stackexchange.com/questions/343308/rasterio-merge-giving-incorrect-result

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Re: Clarification on usage with QGIS

Ratcliff, Christina (A&F, Waite Campus)
 

Ari, I've been running the Geopandas (pandas, shapely fiona), RasterIO configuration since early 2019, and as long as you have access to those gdal backwards compatibility dll's you should be fine.

With libraries depending on GDAL, where possible use those available in osgeo4w. Osgeo4w will then properly install and get them talking to both GDAL and QGIS.

Good luck.


Re: How to specify the name of the input files with the geographical metadata?

juanjo.gomeznavarro@...
 

Thanks. It's a bit strange that this is not a native feature of GDAL, as it seems pretty natural to work on bunches of files that share the same geographical information. But anyway, I have finally opted for using programmatically generated VRT files.


Re: Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()

Patrick Gray
 

Hi Sean and thanks for the reply! I've changed my code to use rasterio.warp.reproject, but the final geotiff is still a rectangle rather than a trapezoid and not actually warping to the GCPs. It is roughly where it should be but not even in the middle of the GCPs which are at each corner of the image. Should I be using that rasterio.transform.from_gcps() function to get the transform for the reproject() function? I must be missing something... Here is my code:

import rasterio
from rasterio.warp import reproject
from rasterio.enums import Resampling
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
        
        print(profile)
        # create the transform from the GCPs
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
        print(np.array(src.read()))
        
        src_crs = "EPSG:4326"  # This is the crs of the GCPs
 
        # Destination: a 1024 x 1024 dataset in WGS84
        dst_raster = np.zeros((1024, 1024), dtype=np.uint16)
        dst_crs = "EPSG:4326"
 
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
 
        new_raster, dst_transform = reproject(
            source=np.array(src.read()),
            destination=dst_raster,
            dst_transform=tsfm,
            gcps=[gcp1,gcp2,gcp3,gcp4],
            src_crs=src_crs,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
            #**kwargs
        )
 
        profile.update(
            dtype=rasterio.uint16,
            transform = dst_transform,
            crs=dst_crs,
            width=dst_raster.shape[0], # TODO unsure if this is correct order
            height=dst_raster.shape[1]
        )
 
        print(profile)
 
        with rasterio.open('example_uint16_tsfm.tif', 'w', **profile) as dst:
            dst.write(src.read(1).astype(rasterio.uint16), 1)


And here is what the image looks like in WGS84:


Re: How to specify the name of the input files with the geographical metadata?

Sean Gillies
 

Hi,

On Tue, Feb 23, 2021 at 9:28 AM <juanjo.gomeznavarro@...> wrote:
Good morning,

I'm working with PNG files, and the metadata is stored in two separate files: .pgw and .png.aux.xml. Currently I'm forced to use the same name for the three files, i.e. the image and the two auxiliary files with the geographical information. This is cumbersome and inconvenient, because I have a set of images that share the same georeferentiation and currently I'm forced to clone the .pgw and .png.aux.xml files in a somewhat silly way for introducing each image. Is there a way to tell rasterio.open that the two files with the metadata are named different from the PNG image?

Thanks a lot!

The naming of the auxiliary files is a GDAL convention and cannot, to my knowledge, be changed. Neither does rasterio allow a user to open an image in read only mode and add or replace the spatial referencing and other metadata. I agree that would be handy in your case. 

Workarounds include:

Opening your PNGs in "r+" mode and patching them, like

with rasterio.open("example.png", "r+") as dataset:
    dataset.crs = rasterio.crs.CRS.from_epsg(3857)
    dataset.transform = rasterio.transform.Affine.translation(ulx, uly) * rasterio.transform.Affine.scale(dx, -dy)

Note that this would produce .pgw and .png.aux.xml files on your system.

Or you could make a generic VRT XML wrapper template and substitute PNG paths into it. See https://gdal.org/drivers/raster/vrt.html. The XML string can then be opened by rasterio.

with rasterio.open("<VRTDataset>...</VRTDataset>") as dataset:
    ...

That's not a well known feature of GDAL/rasterio but I use it quite a bit.

--
Sean Gillies

21 - 40 of 764