Date   
Re: Rasterio 1.0.14

Sean Gillies
 

Yes, you are absolutely right. Google cloud support landed in master and I failed to cherry pick it before 1.0.14 as I'd intended. Sorry! I'm going to release a 1.0.15 to fix this tonight.


On Sun, Jan 27, 2019 at 3:09 PM <madhav@...> wrote:
@Sean
I think something is not right with the PyPI distribution of rasterio 1.0.14. I did pip install rasterio into a clean virtualenv and it installed rasterio. See the pip show rasterio output:

rio
---
Metadata-Version: 2.1
Name: rasterio
Version: 1.0.14
Summary: Fast and direct raster I/O for use with Numpy and SciPy
Author: Sean Gillies
Author-email: sean@...
Installer: pip
License: BSD
Location: /usr/local/lib/python2.7/dist-packages
Requires: cligj, snuggs, click, enum34, attrs, affine, numpy, click-plugins
Classifiers:
  Development Status :: 4 - Beta
  Intended Audience :: Developers
  Intended Audience :: Information Technology
  Intended Audience :: Science/Research
  License :: OSI Approved :: BSD License
  Programming Language :: C
  Programming Language :: Cython
  Programming Language :: Python :: 2.7
  Programming Language :: Python :: 3.3
  Programming Language :: Python :: 3.4
....

However, when I look at site-packages at /usr/local/lib/python2.7/dist-packages/rasterio/ and code for path.py session.py there is no GSSession class or no gs:// paths in the remote schemes in path.py. I got around by git pull source from mapbox/rasterio and doing a local pip install -e . and the above code I pasted works. I don't know if something is wrong with PyPI package distribution. I am curious if others have had success with the latest rasterio 1.0.14 and the new features.
pip show rasterio

pip show rasterio



--
Sean Gillies

Re: Rasterio 1.0.14

Madhav Desetty
 

Great! I am good for now through local rasterio 1.0.14 packaging but happy to test out GS support  again once rasterio 1.0.15 is out.


On Sun, Jan 27, 2019 at 6:15 PM Sean Gillies via Groups.Io <sean=mapbox.com@groups.io> wrote:
Yes, you are absolutely right. Google cloud support landed in master and I failed to cherry pick it before 1.0.14 as I'd intended. Sorry! I'm going to release a 1.0.15 to fix this tonight.

On Sun, Jan 27, 2019 at 3:09 PM <madhav@...> wrote:
@Sean
I think something is not right with the PyPI distribution of rasterio 1.0.14. I did pip install rasterio into a clean virtualenv and it installed rasterio. See the pip show rasterio output:

rio
---
Metadata-Version: 2.1
Name: rasterio
Version: 1.0.14
Summary: Fast and direct raster I/O for use with Numpy and SciPy
Author: Sean Gillies
Author-email: sean@...
Installer: pip
License: BSD
Location: /usr/local/lib/python2.7/dist-packages
Requires: cligj, snuggs, click, enum34, attrs, affine, numpy, click-plugins
Classifiers:
  Development Status :: 4 - Beta
  Intended Audience :: Developers
  Intended Audience :: Information Technology
  Intended Audience :: Science/Research
  License :: OSI Approved :: BSD License
  Programming Language :: C
  Programming Language :: Cython
  Programming Language :: Python :: 2.7
  Programming Language :: Python :: 3.3
  Programming Language :: Python :: 3.4
....

However, when I look at site-packages at /usr/local/lib/python2.7/dist-packages/rasterio/ and code for path.py session.py there is no GSSession class or no gs:// paths in the remote schemes in path.py. I got around by git pull source from mapbox/rasterio and doing a local pip install -e . and the above code I pasted works. I don't know if something is wrong with PyPI package distribution. I am curious if others have had success with the latest rasterio 1.0.14 and the new features.
pip show rasterio

pip show rasterio



--
Sean Gillies



--
Madhav Desetty
Chief Software Architect
Airbus Aerial
(404) 918-0481

Rasterio 1.0.15

Sean Gillies
 

Hi all,

1.0.15 is available on PyPI now and fixes two problems with 1.0.14 and its binary wheels. The first is that 1.0.14 didn't have the advertised Google cloud storage support. The second is that the wheels lacked the advertised HTTP/2 support.

--
Sean Gillies

Re: Rasterio 1.0.15

Madhav Desetty
 

@Sean
Just tested 1.0.15, GS support works great! Thanks...

Madhav

Re: Rasterio 1.0.15

Sean Gillies
 

Hi again,

If you've installed a 1.0.15 rasterio wheel for Linux and are wondering where the HTTP/2 support is, it's in the new 1.0.15-1 wheels on PyPI now.

Contrary to what I announced earlier today, we don't have HTTP/2 support in the OS X wheels. We're relying on the system's libcurl for SSL reasons and it does not have the nghttp2 module.

Thanks for your understanding,

On Mon, Jan 28, 2019 at 9:38 AM Sean Gillies <sean@...> wrote:
Hi all,

1.0.15 is available on PyPI now and fixes two problems with 1.0.14 and its binary wheels. The first is that 1.0.14 didn't have the advertised Google cloud storage support. The second is that the wheels lacked the advertised HTTP/2 support.

--
Sean Gillies


--
Sean Gillies

Re: Rasterio 1.0.15

Madhav Desetty
 

@Sean
I am building a Docker image with the requirements.txt file set to rasterio==1.0.15. Do I have to do anything specific to get the 1.0.15-1 wheels for Docker image running on Linux? 

Re: Rasterio 1.0.15

vincent.sarago@...
 

Hi @Madhav, 
On linux it should work with `rasterio==1.0.15` 

Re: Rasterio 1.0.15

Sean Gillies
 

Correct, there is no need to specify the build tag. There is also no way to override the build tag other than to specify the entire URL of the wheel, but that's a good thing.

On Tue, Jan 29, 2019 at 9:18 AM <vincent.sarago@...> wrote:
Hi @Madhav, 
On linux it should work with `rasterio==1.0.15` 



--
Sean Gillies

New rasterio versions don't handle blank SRS gracefully

nicholas.maxwell@...
 

Short version:

The following causes errors, in rasterio 1.0.15,  but is fine in rasterio 1.0.9. 

from rasterio.crs import CRS
print(CRS())
The following are the only tests in rasterio to check blank CRSs.

def test_null_crs_equality():
assert CRS() == CRS()
a = CRS()
assert a == a
assert not a != a


def test_null_and_valid_crs_equality():
assert (CRS() == CRS(init='epsg:4326')) is False

I feel like print(CRS()) should not cause errors - CRS().__str__ should just return an empty string, for example.

Longer version:

We work with pre-orthorectified sat imagery, which have blank geotransforms (but may have GCPs set), eg,

gdalinfo /data/16APR27183709-P1BS-502376654080_01_P009.tif 
Driver: GTiff/GeoTIFF
Files: /data/16APR27183709-P1BS-502376654080_01_P009.tif
Size is 3201, 3281
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
...

This causes problems on newer rasterio versions (presumably due to the switch to WKT representation for CRS). Example:

import rasterio
with rasterio.open('/data/16APR27183709-P1BS-502376654080_01_P009.tif') as src:
print("CRS:", src.crs)

On rasterio 1.0.9, this is fine, src.crs evaluates to None:

python ~/test.py
/home/nick/venvs/ras9/lib/python3.6/site-packages/rasterio-1.0.9-py3.6-linux-x86_64.egg/rasterio/__init__.py:217: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: None

The NotGeoreferencedWarning is fine, we explicitly ignore it. The point is that src.crs is None.
On 1.0.15, however,

python ~/test.py
/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/__init__.py:216: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: Traceback (most recent call last):
  File "rasterio/_crs.pyx", line 274, in rasterio._crs._CRS.to_dict
  File "rasterio/_err.pyx", line 194, in rasterio._err.exc_wrap_ogrerr
rasterio._err.CPLE_BaseError: OGR Error code 7

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/nick/test.py", line 3, in <module>
    print("CRS:", src.crs)
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 245, in to_string
    return self.to_wkt() or self.to_proj4()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 106, in to_proj4
    return ' '.join(['+{}={}'.format(key, val) for key, val in self.data.items()])
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 170, in data
    self._data = self.to_dict()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 164, in to_dict
    return self._crs.to_dict()
  File "rasterio/_crs.pyx", line 277, in rasterio._crs._CRS.to_dict
rasterio.errors.CRSError: The WKT could not be parsed. OGR Error code 7

Interestingly, the GDAL Python API handles this in a different way:

from osgeo import gdal
ds = gdal.Open('/data/16APR27183709-P1BS-502376654080_01_P009.tif')
prj = ds.GetProjection()
print(prj)

which prints

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
I guess GDAL assumes the GCP projection, in case of blank geotransform.

In any case, I think the src.crs is None behavior is good, since we can just look at src.gcps for the gcps.

I will add a workaround on our side, since
print(CRS().is_valid)
works, (that is, CRS().is_valid evaluates to False, without throwing exceptions.

Re: New rasterio versions don't handle blank SRS gracefully

Sean Gillies
 

Hi Nicholas,

Thanks for the explanation. We have at least one more related regression in 1.0.15 and I intend to fix both of these in 1.0.16.


On Fri, Feb 1, 2019 at 10:36 AM <nicholas.maxwell@...> wrote:
Short version:

The following causes errors, in rasterio 1.0.15,  but is fine in rasterio 1.0.9. 

from rasterio.crs import CRS
print(CRS())
The following are the only tests in rasterio to check blank CRSs.

def test_null_crs_equality():
assert CRS() == CRS()
a = CRS()
assert a == a
assert not a != a


def test_null_and_valid_crs_equality():
assert (CRS() == CRS(init='epsg:4326')) is False

I feel like print(CRS()) should not cause errors - CRS().__str__ should just return an empty string, for example.

Longer version:

We work with pre-orthorectified sat imagery, which have blank geotransforms (but may have GCPs set), eg,

gdalinfo /data/16APR27183709-P1BS-502376654080_01_P009.tif 
Driver: GTiff/GeoTIFF
Files: /data/16APR27183709-P1BS-502376654080_01_P009.tif
Size is 3201, 3281
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
...

This causes problems on newer rasterio versions (presumably due to the switch to WKT representation for CRS). Example:

import rasterio
with rasterio.open('/data/16APR27183709-P1BS-502376654080_01_P009.tif') as src:
print("CRS:", src.crs)

On rasterio 1.0.9, this is fine, src.crs evaluates to None:

python ~/test.py
/home/nick/venvs/ras9/lib/python3.6/site-packages/rasterio-1.0.9-py3.6-linux-x86_64.egg/rasterio/__init__.py:217: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: None

The NotGeoreferencedWarning is fine, we explicitly ignore it. The point is that src.crs is None.
On 1.0.15, however,

python ~/test.py
/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/__init__.py:216: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: Traceback (most recent call last):
  File "rasterio/_crs.pyx", line 274, in rasterio._crs._CRS.to_dict
  File "rasterio/_err.pyx", line 194, in rasterio._err.exc_wrap_ogrerr
rasterio._err.CPLE_BaseError: OGR Error code 7

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/nick/test.py", line 3, in <module>
    print("CRS:", src.crs)
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 245, in to_string
    return self.to_wkt() or self.to_proj4()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 106, in to_proj4
    return ' '.join(['+{}={}'.format(key, val) for key, val in self.data.items()])
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 170, in data
    self._data = self.to_dict()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/crs.py", line 164, in to_dict
    return self._crs.to_dict()
  File "rasterio/_crs.pyx", line 277, in rasterio._crs._CRS.to_dict
rasterio.errors.CRSError: The WKT could not be parsed. OGR Error code 7

Interestingly, the GDAL Python API handles this in a different way:

from osgeo import gdal
ds = gdal.Open('/data/16APR27183709-P1BS-502376654080_01_P009.tif')
prj = ds.GetProjection()
print(prj)

which prints

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
I guess GDAL assumes the GCP projection, in case of blank geotransform.

In any case, I think the src.crs is None behavior is good, since we can just look at src.gcps for the gcps.

I will add a workaround on our side, since
print(CRS().is_valid)
works, (that is, CRS().is_valid evaluates to False, without throwing exceptions.



--
Sean Gillies

Re: Implementing -scale functionality of gdal_translate

linden.kyle@...
 

Hi Sean, I'm having some trouble with the implementation of this and was wondering if you might be able to give me a hand. Here's the code I am working with: https://pastebin.com/kecx0D2b. What I'm trying to do is take a geotiff and convert it to a set of TMS-formatted tiles (including reprojecting to web mercator). In this example `tile_bounds` are the bounds of the tile in web mercator. I originally had this code sample working until I attempted to add in the code to do the pixel scaling. Here's the original working code: https://pastebin.com/dZVRVKi2

Here's what the output looks like with the broken code: https://imgur.com/xLtWrGZ
And here's the correct version (but no support for pixel scaling): https://imgur.com/aoNCyIp

Reading an image by overlapped blocks of chosen depth

mail2abhisek.maiti@...
 

In rasterio, reading image array in the unit of blocks is pretty convenient and useful. However, Overlapping blocks could be very useful for certain applications which requires kernel wise processing (such as simple smoothing filter).  So far there is no provision to make a block overlap with its neighboring blocks as per requirement of the application. It would be immensely handy if this feature gets implemented in rasterio. Related Stackoverflow Question.

Re: Implementing -scale functionality of gdal_translate

linden.kyle@...
 

Ok, so it appears the issue was that I was attempting to do a windowed read outside the bounds of the geotiff (for TMS edge tiles) and `src.read()` was only returning the data inside the bounds since I wasn't using the `boundless=True` option. However, the next problem I am running into is that I'd like to use `fill_value=numpy.nan` to fill in the edge data with NaN values. However, my data is UInt8. Is there any way to specify an option to cast the data read by the call to `read()` into float32 so that I can support this?

Re: Implementing -scale functionality of gdal_translate

Sean Gillies
 

Hi,

I apologize for being slow to reply, I've been occupied with fixing bugs in the last couple of releases.

In theory it is possible to add an option to the read method which will cast the data because GDAL already supports this well in versions 2 and up. I haven't pursued this because converting unsigned int to float is trivial, if not cheap, with numpy:

   data = dataset.read().astype('float64')

And if you pass masked=True to read, you'll get a masked array and then can apply a fill_value of numpy.nan under the invalid data mask. Is this useful for you or have you already been down this path?

GDAL's data type casting would not produce an extra copy of the data and so would be optimal for some applications.


On Tue, Feb 5, 2019 at 11:09 AM <linden.kyle@...> wrote:
Ok, so it appears the issue was that I was attempting to do a windowed read outside the bounds of the geotiff (for TMS edge tiles) and `src.read()` was only returning the data inside the bounds since I wasn't using the `boundless=True` option. However, the next problem I am running into is that I'd like to use `fill_value=numpy.nan` to fill in the edge data with NaN values. However, my data is UInt8. Is there any way to specify an option to cast the data read by the call to `read()` into float32 so that I can support this?



--
Sean Gillies

Re: Implementing -scale functionality of gdal_translate

linden.kyle@...
 

Sean, no problem. Thanks for writing back. I'll keep that in mind for the future. What I ended up doing instead was modifying the tile window transform (which is passed to the reproject method) in the case where either the left or top coordinate is outside the bounds of the geotiff. In that case I change the translation of the transform to the left and/or top bound of the geotiff respectively. Not sure if this is the "right" way to do it or not, but it seems to be working.

Re: Implementing -scale functionality of gdal_translate

linden.kyle@...
 

So I'm seeing some odd behavior now. I ended up trying the masked_array solution, but I'm seeing breaks in my data at the tile boundaries.

Here's what the data looks like if I use reprojection with just the source band (code: https://pastebin.com/dZVRVKi2): https://imgur.com/WpFKfij

Here's what the data looks like I read into a masked array and use an ndarray source (code: https://pastebin.com/z1xd8mXu): https://imgur.com/2QkQQZf
Ironically enough I get the exact same result with my solution that doesn't use a masked array, but instead modifies the transform: https://pastebin.com/7ph1VYWX

Re: Reading an image by overlapped blocks of chosen depth

Luke Pinner
 

I've put a couple of methods in my answer to your SO question, by expanding the block windows and also for arbitrary windows (unrelated to internal blocks/tiles).

Re: Reading an image by overlapped blocks of chosen depth

mail2abhisek.maiti@...
 

I was not aware of the boundless reading. Fantastic implementation. Thanks.

Rasterio 1.0.18

Sean Gillies
 

Hi all,

We've made three new releases this week to fix regressions introduced in version 1.0.14. Please upgrade to 1.0.18 to take advantage of the changes in 1.0.14, but without suffering from the bugs. See https://github.com/mapbox/rasterio/blob/master/CHANGES.txt for more details.

--
Sean Gillies

Re: Implementing -scale functionality of gdal_translate

Sean Gillies
 

Hi,

On Wed, Feb 6, 2019 at 2:32 PM <linden.kyle@...> wrote:
So I'm seeing some odd behavior now. I ended up trying the masked_array solution, but I'm seeing breaks in my data at the tile boundaries.

Here's what the data looks like if I use reprojection with just the source band (code: https://pastebin.com/dZVRVKi2): https://imgur.com/WpFKfij

Here's what the data looks like I read into a masked array and use an ndarray source (code: https://pastebin.com/z1xd8mXu): https://imgur.com/2QkQQZf
Ironically enough I get the exact same result with my solution that doesn't use a masked array, but instead modifies the transform: https://pastebin.com/7ph1VYWX

Your code looks okay. I recommend padding your read window by a few pixels on each side so that tile_window, tile_data, and tile_transform are slight dilated. I think this is very likely to eliminate what looks like an off-by-one error in the reprojection, an error that you don't see when you use the entire source dataset as the reprojection source. Does that make sense?

--
Sean Gillies