Date   

Re: Should rasterio be compatible with numpy.ma?

Sean Gillies
 

Hi,

I'm in favor of making filling the default for a masked array. In the same vein, I wonder if write should also call write_mask if a masked array is passed?


On Thu, Nov 11, 2021 at 8:54 AM <dreaminghexnut@...> wrote:
I was just tripped up by the fact that rasterio (I'm on 1.2.10) ignores the mask in numpy masked arrays. It simply writes the masked array's data array. If that underlying data array has the same fill value as the one defined in the profile, that works as expected. However, if the fill value in the data array is different from the one in the profile or if there are original brightness values instead of a fill value, the mask is completely ignored.

There's a simple workaround:
```
meta.update(fill_value=masked_img.fill_value)
with open(target, "w", **meta) as dst:
    dst.write(masked_img.filled())
```

But shouldn't that be the default behaviour? That rasterio should call .filled() every time it receives a Numpy masked array?

Cheers,
JOey



--
Sean Gillies


Should rasterio be compatible with numpy.ma?

dreaminghexnut@...
 

I was just tripped up by the fact that rasterio (I'm on 1.2.10) ignores the mask in numpy masked arrays. It simply writes the masked array's data array. If that underlying data array has the same fill value as the one defined in the profile, that works as expected. However, if the fill value in the data array is different from the one in the profile or if there are original brightness values instead of a fill value, the mask is completely ignored.

There's a simple workaround:
```
meta.update(fill_value=masked_img.fill_value)
with open(target, "w", **meta) as dst:
    dst.write(masked_img.filled())
```

But shouldn't that be the default behaviour? That rasterio should call .filled() every time it receives a Numpy masked array?

Cheers,
JOey


Re: extending rio-calc

Gregory, Matthew
 

Thanks Sean,

 

This is really helpful.  It works for the handful of functions in my workflow, so this will get me through.  It’s always humbling to find out a corner of Python that I didn’t know existed - thanks for pointing me to sitecustomize.py.

 

Not that I would assume that you’d want to manage this, but is there a pathway (or desire) for user-contributed functions to make their way into snuggs/rio-calc or would you prefer the user customization route?

 

thanks, matt

 

 

From: main@rasterio.groups.io <main@rasterio.groups.io> On Behalf Of Sean Gillies via groups.io
Sent: Saturday, October 30, 2021 11:40 AM
To: main@rasterio.groups.io
Subject: Re: [rasterio] extending rio-calc

 

[This email originated from outside of OSU. Use caution with links and attachments.]

Hi Matt,

 

No, we don't have a plugin system for rio-calc. If this is for your personal workflows, maybe a sitecustomize.py or usercustomize.py file that imported and patched snuggs would be worth a try?

 

I saved the following to sitecustomize.py in the site-packages directory of a Python environment I named "play".

 

import snuggs
snuggs.func_map["test"] = lambda *args: print(*args)

 

This isn't going to really work because it doesn't return any values that rio-calc can use, but you can see that the module does get imported and snuggs does get patched:

 

(play) $ rio calc "(test 'hi hi')" RGB.byte.tif -o /tmp/calc.tif
hi hi
Traceback (most recent call last):

   ...

 

This is the closest thing Python has to built-in extensibility. Limitations are discussed at https://www.python.org/dev/peps/pep-0648/#limitations-of-sitecustomize-py.

 

 

On Fri, Oct 29, 2021 at 8:45 AM Gregory, Matthew <matt.gregory@...> wrote:

Hi all,

I've found rio-calc to be really useful for raster algebra expressions.  I've come up with a few additional functions that are useful in my workflow that seem to work well within rio-calc.

I'm certain, however, that I don't know how to introduce these functions into rio-calc without directly registering these functions with snuggs within calc.py (e.g. https://github.com/rasterio/rasterio/blob/master/rasterio/rio/calc.py#L159-L164). 

Is there a plug-in way of extending rio-calc to bring in additional user-defined functions?

thanks, matt






--

Sean Gillies


Re: extending rio-calc

Sean Gillies
 

Hi Matt,

No, we don't have a plugin system for rio-calc. If this is for your personal workflows, maybe a sitecustomize.py or usercustomize.py file that imported and patched snuggs would be worth a try?

I saved the following to sitecustomize.py in the site-packages directory of a Python environment I named "play".

import snuggs
snuggs.func_map["test"] = lambda *args: print(*args)

This isn't going to really work because it doesn't return any values that rio-calc can use, but you can see that the module does get imported and snuggs does get patched:

(play) $ rio calc "(test 'hi hi')" RGB.byte.tif -o /tmp/calc.tif
hi hi
Traceback (most recent call last):
   ...

This is the closest thing Python has to built-in extensibility. Limitations are discussed at https://www.python.org/dev/peps/pep-0648/#limitations-of-sitecustomize-py.


On Fri, Oct 29, 2021 at 8:45 AM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

I've found rio-calc to be really useful for raster algebra expressions.  I've come up with a few additional functions that are useful in my workflow that seem to work well within rio-calc.

I'm certain, however, that I don't know how to introduce these functions into rio-calc without directly registering these functions with snuggs within calc.py (e.g. https://github.com/rasterio/rasterio/blob/master/rasterio/rio/calc.py#L159-L164). 

Is there a plug-in way of extending rio-calc to bring in additional user-defined functions?

thanks, matt







--
Sean Gillies


extending rio-calc

Gregory, Matthew
 

Hi all,

I've found rio-calc to be really useful for raster algebra expressions. I've come up with a few additional functions that are useful in my workflow that seem to work well within rio-calc.

I'm certain, however, that I don't know how to introduce these functions into rio-calc without directly registering these functions with snuggs within calc.py (e.g. https://github.com/rasterio/rasterio/blob/master/rasterio/rio/calc.py#L159-L164).

Is there a plug-in way of extending rio-calc to bring in additional user-defined functions?

thanks, matt


Re: mapbox/rasterio -> rasterio/rasterio

Sean Gillies
 

The transfer is complete.

On Thu, Oct 21, 2021 at 8:51 AM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
Heads up. everybody,

I've been wanting to move rasterio to its own GitHub organization for some time and it is finally going to happen. Last week the Mapbox legal team signed off on transferring rasterio to its own GitHub org and today the IT team is going to transfer it to https://github.com/rasterio. I apologize for the short notice. You may want to pause development until the transfer is complete.

I expect the transfer to be nearly seamless. For once, the Python Package Index's lack of project namespaces is a feature, not a bug. As a user, you won't have to change any of your dependency specifications :)


--
Sean Gillies


mapbox/rasterio -> rasterio/rasterio

Sean Gillies
 

Heads up. everybody,

I've been wanting to move rasterio to its own GitHub organization for some time and it is finally going to happen. Last week the Mapbox legal team signed off on transferring rasterio to its own GitHub org and today the IT team is going to transfer it to https://github.com/rasterio. I apologize for the short notice. You may want to pause development until the transfer is complete.

I expect the transfer to be nearly seamless. For once, the Python Package Index's lack of project namespaces is a feature, not a bug. As a user, you won't have to change any of your dependency specifications :)
 
--
Sean Gillies


Rasterio 1.3a2

Sean Gillies
 

Hi all,

The second 1.3 pre-release is on the package index this morning. There's one new feature: PROJJSON support (thank you, Alan Snow) and one change to the reprojection code that prevents unnecessary data copying when warping raster data stored in a numpy array (thank you, Ryan Grout). Your feedback on these would be greatly appreciated!


Re: Warp doesn't work with a general oblique transform?

erik.koene@...
 

Hi Denis, I'm ashamed and happy to say that that indeed did the trick, the warp works perfectly now. Thanks!


Re: concurrent.futures errors running rio mbtiles on Windows

Sean Gillies
 

Hi Matt,

I'm sorry that you're having trouble. We don't test rigorously on Windows and so it's not surprising that some aspects of rasterio aren't correct.

If we execute the tiler as `rio mbtiles`, the __main__ module is generated by setuptools and my understanding is that this hides, on Windows, imports needed by workers.

Could you try adding, after https://github.com/mapbox/rio-mbtiles/blob/main/mbtiles/scripts/cli.py#L600, the following block?

if __name__ == "__main__":
    mbtiles()

And then try executing the program not as a rio plugin, but as python -m mbtiles.scripts.cli ... ?


On Mon, Oct 11, 2021 at 11:36 AM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

Preface this by saying I know next to nothing about multiprocessing and my debugging skills are failing me.

I'm trying to run a simple rio mbtiles command on Windows (Python 3.8, rasterio 1.2.9, rio_mbtiles 1.6.0) with a small test GeoTiff in Web Mercator using this syntax:

  > rio mbtiles
      --resampling average
      --zoom-levels 5..8
      --format PNG
      -j 4
      --progress-bar
      test.tif
      test.mbtiles

I often get this error stack trace:

  """
  File "c:\users\gregorma\envs\default-38\lib\site-packages\mbtiles\cf.py",
  line 62, in process_tiles
    tile, contents = future.result()
  File "c:\users\gregorma\appdata\local\programs\python\python38\
  lib\concurrent\futures\_base.py", line 432, in result
    return self.__get_result()
  File "c:\users\gregorma\appdata\local\programs\python\python38\
  lib\concurrent\futures\_base.py", line 388, in __get_result
    raise self._exception
  concurrent.futures.process.BrokenProcessPool: A process in the process pool
  was terminated abruptly while the future was running or pending.
  """

I can set up the same environment using WSL (Debian) and run the same command and it works flawlessly.  I don't think it's relevant but I'm using virtualenvs on both Windows and WSL.

I've read numerous posts about including either the if __name__ == "__main__" guard or (with multiprocessing) using freeze_support for helping with these issues on Windows, but have no idea where these would be inserted for this command (either at rio itself or within any of the plugins that use concurrent.futures). 

Has anyone had luck running rio mbtiles on a Windows box?  I obviously have the WSL workaround if needed.  Please let me know if more information is needed.

Thanks for any advice,
matt









--
Sean Gillies


Re: Warp doesn't work with a general oblique transform?

Denis Rykov
 

I would recommend you to upgrade rasterio. The version you are using is quite old.


Re: Why does Dataset.transform return the identity matrix if there is no geotransform?

Sean Gillies
 

Hi,

Yes, we're following GDAL, which returns an identity matrix instead of NULL. https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd.

GDAL's design decision has complicated things for rasterio. See the discussion in https://github.com/mapbox/rasterio/issues/524 for example.


On Mon, Oct 18, 2021 at 8:18 AM <interdimensional.cabbage@...> wrote:
I have been reading -> modifying -> writing some GeoTiffs and I got bit by Dataset.transform returning the identity transform. There's a big difference between there being no geotransform and there being the identity geotransform.

The result of this discrepancy is that if if you simply try to copy the transform from the source dataset to the destination dataset, it now has completely incorrect information. There are work arounds; you can detect that it's the identity and set it back to None. But this seems like an unnecessary asymmetry in reading/writing.

I was thinking of raising an issue on GitHub, but I figured I'd ask first. Is there any reason that Dataset.transform should return Affine.identity() instead of None? Is it to match GDAL in some way?



--
Sean Gillies


Re: Warp doesn't work with a general oblique transform?

erik.koene@...
 

Hi Sean, thanks for your response. I'm running Python 3.7, rasterio version 1.0.14; gdalinfo reports version 2.4, and proj reports version 5.2.0. I'd like to clarify with you that this is a result from rasterio.crs.CRS.from_string, and not from another library? Because when I use pyproj, the CRS is in principle created correctly,

So, question 1: how do you generate the CRS, with pyproj or with rasterio? Below, you can see that the pyproj call returns no problem yet, while the rasterio version already does.
import pyproj
dst_crs = pyproj.CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')
print(dst_crs) # >> +proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137 +type=crs
import rasterio
dst_crs = rasterio.crs.CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')

Question 2, when I use rasterio.crs.CRS.from_wkt() and use the string you post, then the warp still returns an error:
---------------------------------------------------------------------------
CRSError                                  Traceback (most recent call last)
<ipython-input-6-8658800cad11> in <module>
     29 
     30 with rasterio.open(path) as src:
---> 31     with WarpedVRT(src, **vrt_options) as vrt:
     32         # Read all data into memory and show plot
     33         data = vrt.read()

rasterio/_warp.pyx in rasterio._warp.WarpedVRTReaderBase.__init__()

rasterio/_base.pyx in rasterio._base._osr_from_crs()

CRSError: Invalid CRS: CRS.from_wkt('GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["unknown",METHOD["PROJ ob_tran o_proj=longlat"],PARAMETER["lon_0",-170,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lon_p",-180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lat_p",137,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS["ellipso
It complains about an "invalid" CRS at the stage of the virtual warp, even though the WKT string is accepted.

The oblique transform works fine to transform individual coordinates. I.e., running the following snippet
from pyproj import Transformer
rotated_to_CRS = Transformer.from_crs('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137',
'EPSG:4326',
always_xy=True)
coordinates=[(0,0), (1,1)]
[pt for pt in rotated_to_CRS.itransform(coordinates)] # >> [(9.99999999999999, 47.0), (11.49406628349866, 47.990464430559726)]
correctly transforms from the oblique rotated coordinate system into the WGS84 system. So I don't understand what is strictly wrong with the CRS.


Why does Dataset.transform return the identity matrix if there is no geotransform?

Brandon Victor
 

I have been reading -> modifying -> writing some GeoTiffs and I got bit by Dataset.transform returning the identity transform. There's a big difference between there being no geotransform and there being the identity geotransform.

The result of this discrepancy is that if if you simply try to copy the transform from the source dataset to the destination dataset, it now has completely incorrect information. There are work arounds; you can detect that it's the identity and set it back to None. But this seems like an unnecessary asymmetry in reading/writing.

I was thinking of raising an issue on GitHub, but I figured I'd ask first. Is there any reason that Dataset.transform should return Affine.identity() instead of None? Is it to match GDAL in some way?


Re: Warp doesn't work with a general oblique transform?

Sean Gillies
 

Hi,

With the latest release of GDAL and PROJ 8 I see

>>> CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')
CRS.from_wkt('GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["unknown",METHOD["PROJ ob_tran o_proj=longlat"],PARAMETER["lon_0",-170,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lon_p",-180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lat_p",137,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]')

What versions of rasterio, GDAL, and PROJ do you have? And where did you install them from?


On Wed, Oct 13, 2021 at 10:04 AM <erik.koene@...> wrote:
The following Python code...

import rasterio
from rasterio.crs import CRS
from rasterio import shutil as rio_shutil
from rasterio.vrt import WarpedVRT
import matplotlib.pyplot as plt

path = 'CORINE/U2018_CLC2018_V2020_20u1.tif'

# Destination CRS is 4326
# dst_crs = CRS.from_epsg(4326) # <--- This one will work...
dst_crs = CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137') # <--- This one doesn't work

# Output image dimensions
dst_height = dst_width = 1000

# Output image transform
left, bottom, right, top = 10.897971,49.059812, 10.907593,49.064142 # Coordinates in target CRS
xres = (right - left) / dst_width
yres = (top - bottom) / dst_height
dst_transform = rasterio.Affine(xres, 0.0, left,
                              0.0, -yres, top)

vrt_options = {
    'crs': dst_crs,
    'transform': dst_transform,
    'height': dst_height,
    'width': dst_width,
}

with rasterio.open(path) as src:
    with WarpedVRT(src, **vrt_options) as vrt:
        # Read all data into memory and show plot
        data = vrt.read()
        plt.imshow(data.squeeze())


fails to run, giving the error message
The PROJ4 dict could not be understood. OGR Error code 5
which doesn't appear if I use a more 'standard' destination CRS, e.g., CRS.from_epsg(4326) . I'll welcome any suggestions! Thanks.



--
Sean Gillies


Rasterio 1.3a1

Sean Gillies
 

Hi all,

We have a first alpha pre-release on PyPI now. Many of the changes are noted at https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L7. It's not quite complete, so please see https://github.com/mapbox/rasterio/milestone/103 for what new features are or are not yet in this release. Please remember that some things may be removed or altered a bit before 1.3.0. I'm eager to hear whether rasterio's merge tool is working better in 1.3a1.

--
Sean Gillies


Warp doesn't work with a general oblique transform?

erik.koene@...
 

The following Python code...

import rasterio
from rasterio.crs import CRS
from rasterio import shutil as rio_shutil
from rasterio.vrt import WarpedVRT
import matplotlib.pyplot as plt

path = 'CORINE/U2018_CLC2018_V2020_20u1.tif'

# Destination CRS is 4326
# dst_crs = CRS.from_epsg(4326) # <--- This one will work...
dst_crs = CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137') # <--- This one doesn't work

# Output image dimensions
dst_height = dst_width = 1000

# Output image transform
left, bottom, right, top = 10.897971,49.059812, 10.907593,49.064142 # Coordinates in target CRS
xres = (right - left) / dst_width
yres = (top - bottom) / dst_height
dst_transform = rasterio.Affine(xres, 0.0, left,
                              0.0, -yres, top)

vrt_options = {
    'crs': dst_crs,
    'transform': dst_transform,
    'height': dst_height,
    'width': dst_width,
}

with rasterio.open(path) as src:
    with WarpedVRT(src, **vrt_options) as vrt:
        # Read all data into memory and show plot
        data = vrt.read()
        plt.imshow(data.squeeze())


fails to run, giving the error message
The PROJ4 dict could not be understood. OGR Error code 5
which doesn't appear if I use a more 'standard' destination CRS, e.g., CRS.from_epsg(4326) . I'll welcome any suggestions! Thanks.


Rasterio 1.2.10

Sean Gillies
 

Hi all,

Rasterio has just one bug fix, but it's an important one: https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L15-L19.

The libraries in the wheels are the same as in 1.2.9. We don't have wheels for Python 3.10 yet, but maybe later in the week. Availability of 3.10 on CI platforms is uneven.

--
Sean Gillies


concurrent.futures errors running rio mbtiles on Windows

Gregory, Matthew
 

Hi all,

Preface this by saying I know next to nothing about multiprocessing and my debugging skills are failing me.

I'm trying to run a simple rio mbtiles command on Windows (Python 3.8, rasterio 1.2.9, rio_mbtiles 1.6.0) with a small test GeoTiff in Web Mercator using this syntax:

> rio mbtiles
--resampling average
--zoom-levels 5..8
--format PNG
-j 4
--progress-bar
test.tif
test.mbtiles

I often get this error stack trace:

"""
File "c:\users\gregorma\envs\default-38\lib\site-packages\mbtiles\cf.py",
line 62, in process_tiles
tile, contents = future.result()
File "c:\users\gregorma\appdata\local\programs\python\python38\
lib\concurrent\futures\_base.py", line 432, in result
return self.__get_result()
File "c:\users\gregorma\appdata\local\programs\python\python38\
lib\concurrent\futures\_base.py", line 388, in __get_result
raise self._exception
concurrent.futures.process.BrokenProcessPool: A process in the process pool
was terminated abruptly while the future was running or pending.
"""

I can set up the same environment using WSL (Debian) and run the same command and it works flawlessly. I don't think it's relevant but I'm using virtualenvs on both Windows and WSL.

I've read numerous posts about including either the if __name__ == "__main__" guard or (with multiprocessing) using freeze_support for helping with these issues on Windows, but have no idea where these would be inserted for this command (either at rio itself or within any of the plugins that use concurrent.futures).

Has anyone had luck running rio mbtiles on a Windows box? I obviously have the WSL workaround if needed. Please let me know if more information is needed.

Thanks for any advice,
matt


Re: How can I set config options when using the command-line tools?

robin@...
 

Hi Sean,

That seems to work now - I think I had some problems with how I was setting my environment variables on Windows. Setting them through the GUI configuration dialog (under System->Environment Variables) seems to have worked.

Thanks,

Robin

61 - 80 of 945