satellogic / telluric Goto Github PK

View Code? Open in Web Editor NEW

87.0 14.0 18.0 44.36 MB

telluric is a Python library to manage vector and raster geospatial data in an interactive and easy way

License: MIT License

Makefile 0.01% Jupyter Notebook 89.59% Python 10.37% Dockerfile 0.03%

python-library python python3 maps geospatial-data raster

telluric's Introduction

telluric

Overview

telluric is a Python library to manage vector and raster geospatial data in an interactive and easy way.

Opening a raster is as simple as:

In [1]: import telluric as tl

In [2]: tl.GeoRaster2.open("http://download.osgeo.org/geotiff/samples/usgs/f41078e1.tif")
Out[2]: <telluric.georaster.GeoRaster2 at 0x7facd183ad68>

And reading some vector data is equally simple:

In [3]: tl.FileCollection.open("shapefiles/usa-major-cities.shp")
Out[3]: <telluric.collections.FileCollection at 0x7facd1183048>

For more usage examples and a complete API reference, check out our documentation on Read the Docs.

The source code and issue tracker are hosted on GitHub, and all contributions and feedback are more than welcome.

Installation

You can install telluric using pip:

pip install telluric[vis]

Read more complete installation instructions at our documentation.

telluric is a pure Python library, and therefore should work on Linux, OS X and Windows provided that you can install its dependencies. If you find any problem, please open an issue and we will take care of it.

Development

telluric is usually developed on Linux. For full tests do:

$ make build
$ make test

for testing single tests do:

$ make dockershell
docker$ python -m pytest TEST_FILE::TEST_NAME

Support

Join our Matrix chat to ask all sorts of questions!

telluric's People

Contributors

Stargazers

Watchers

Forkers

astrojuanlu gadgetsteve drnextgis luciotorre amitaronovitch jakop345 nojuman adivicco47 fisadev kirosg florianpignol maralonso eecarres lilsnore enomis-dev jmarintur agumartina pepgonfaus

telluric's Issues

GeoRaster2 `equal` is not always right

There are cases when two rasters have the same information, but comparing them returns False, for example:

two rasters with the same information and different CRS
two rasters with opposite affine scale and transformed image data

Here is a code example:

import telluric as tl
import numpy as np
from affine import Affine
from telluric.constants import WEB_MERCATOR_CRS
from shapely.geometry import Point, Polygon

array = np.arange(0, 20).reshape(1, 4, 5)
array2 = np.arange(19, -1, -1).reshape(1, 4, 5)
array2.sort()

image1 = np.ma.array(array, mask=False)
image2 = np.ma.array(array2, mask=False)

aff2 = Affine.translation(0, -8) * Affine.scale(2, 2)
aff = Affine.scale(2, -2)

r1 = tl.GeoRaster2(image=image1, affine=aff, crs=WEB_MERCATOR_CRS)
r2 = tl.GeoRaster2(image=image2, affine=aff2, crs=WEB_MERCATOR_CRS)

r1 == r2  # should be true in my opinion
r1.to_png() == r2.to_png()  # is True

# comparing different points of the rasters shows they are the same
np.array_equal(r1.get(tl.GeoVector(Point(2, -2), crs=WEB_MERCATOR_CRS)),
               r2.get(tl.GeoVector(Point(2, -2), crs=WEB_MERCATOR_CRS)))   # is True

np.array_equal(r1.get(tl.GeoVector(Point(0, -0), crs=WEB_MERCATOR_CRS)),
               r2.get(tl.GeoVector(Point(0, -0), crs=WEB_MERCATOR_CRS)))  # is True

np.array_equal(r1.get(tl.GeoVector(Point(9, -2), crs=WEB_MERCATOR_CRS)), 
              r2.get(tl.GeoVector(Point(9, -2), crs=WEB_MERCATOR_CRS)))  # is True

np.array_equal(r1.get(tl.GeoVector(Point(9, -7), crs=WEB_MERCATOR_CRS)), 
              r2.get(tl.GeoVector(Point(9, -7), crs=WEB_MERCATOR_CRS)))  # is True

roi = tl.GeoVector(Polygon.from_bounds(0, 0, 3, -3), crs=WEB_MERCATOR_CRS)

r1c = r1.crop(roi)
r2c = r2.crop(roi)

r1c == r2c  # should be true
r2c.to_png() == r1c.to_png()  #  is True

GeoRaster2.align_raster_to_mercator_tiles is broken

In [1]: import telluric as tl

In [2]: r = tl.GeoRaster2.open('tests/data/raster/rgb.tif')

In [3]: r.align_raster_to_mercator_tiles()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-2b64356334c5> in <module>()
----> 1 r.align_raster_to_mercator_tiles()

~/sandbox/telluric/telluric/telluric/georaster.py in align_raster_to_mercator_tiles(self)
   1462         # this requires geographical crs
   1463         gv = raster.footprint().reproject(DEFAULT_CRS)
-> 1464         bounding_tile = mercantile.bounding_tile(*gv.shape.bounds)
   1465         window = raster._tile_to_window(*bounding_tile)
   1466         width = math.ceil(abs(window.width))

~/sandbox/telluric/telluric/telluric/vectors.py in __getattr__(self, item)
    251 
    252         else:
--> 253             raise AttributeError("'{}' object has no attribute '{}'".format(self.__class__, item))
    254 
    255         # Return the newly bound attribute

AttributeError: '<class 'telluric.vectors.GeoVector'>' object has no attribute 'shape'

Also it is not clear where is _tile_to_window method comes from.

Get_tile should return a raster which is projected to web_mercator

Add conda-forge packages

Documentation: https://conda-forge.org/docs/recipe.html

Documentation: Area is in square meters not meters

On the readthedocs page: http://telluric.readthedocs.io/en/v0.1.0/User%20Guide.html

Expected behavior and actual behavior

The comment on In [6]: gv1.area # Real area in meters this should read In [6]: gv1.area # Real area in square meters

Steps to reproduce the problem

Just read the documents.

Operating system and installation details

Any & All.

Collaboration with PyTroll

Hi,

I couldn't find any information in your documentation about a chat room or mailing list so I thought I'd make an issue here. I'm a core developer on the PyTroll team (http://pytroll.github.io/). I saw your email on the python-announce mailing list and thought we could discuss collaboration. Your package is very interesting and you have some similar features to our package SatPy (https://github.com/pytroll/satpy) but you seem to lean on the GIS side of things where SatPy is currently focused on reading satellite raster data from various satellite imager data files. You may also be interested in our PyResample library (https://github.com/pytroll/pyresample) which can help resample data to uniform grids (a.k.a. areas). Lastly your orbit-predictor library seems very similar to our pyorbital package (https://github.com/pytroll/pyorbital).

I'm not sure where or how we can collaborate right now, but wanted to make our projects known to each other so we can collaborate in the future if/when it makes sense. I'll also tag @mraspaud here who is another core developer on the PyTroll team.

get_tile - float division by zero

import telluric as tl
url = "https://tellurictest.blob.core.windows.net/test/dataset240/2538/bad.tif?sp=rl&st=2018-07-09T06:15:52Z&se=2018-07-10T06:15:52Z&sv=2017-11-09&sig=Vjkihrlxqn84F2UwFBQc8Sn0UEgyzZ2%2FPJOh%2Bl6HP%2FM%3D&sr=b"
r = tl.GeoRaster2.open(url)
r.get_tile(2018,1567,12)

outputs:

ZeroDivisionError: float division by zero

Telluric versions:
latest, master

factors argument is not respected

factors argument of GeoRaster2.save is not respected after 2ec9194. Also the part of this doc related to the default value is not actual any more:

telluric/telluric/georaster.py

Line 535 in d7206c3

    
                   * factors: list of factors for the overview, default: [2, 4, 8, 16, 32, 64, 128]

Remove nodata_value from rasterize

We should remove the nodata_value parameter from rasterization functions: telluric always works with masked arrays anyway, so this should remain as an implementation detail. We should still do the right thing in all cases (even when the user specifies fill_value=0).

cc @luciotorre

Write down development process

Ideas after telluric 0.2.0 is released (#102):

Automate release notes using https://github.com/hawkowl/towncrier (see https://pip.pypa.io/en/stable/development/#adding-a-news-entry)
Explicitly use trunk-based development (https://www.toptal.com/software/trunk-based-development-git-flow)
- IMHO, the first two criticism points of trunk-based do not really apply because changes must be approved by maintainers
Explain scope, vision, core dependencies, roadmap
Write down governance model

Debug mode

We should have debug mode that will help us understand what is going on.
In this mode we should see:

How much time take computation
When using 3rd parties what are the arguments that are passed to them
...

slicing a FeatureCollection returns a list

Expected behavior and actual behavior

Slicing a FeatureCollection should return another FeatureCollection, but now it returns a list.

Steps to reproduce the problem

In [4]: fc = tl.FeatureCollection([
   ...: tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr1': 1}),
   ...: tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr2': '2'})
   ...: ])

In [5]: fc
Out[5]: <telluric.collections.FeatureCollection at 0x7f9f949d0908>

In [6]: fc[:1]
Out[6]: [GeoFeature(Point, {'attr1': 1})]

Operating system and installation details

Linux, latest telluric.

Reduce default zoom for points

When one displays a point in telluric, the zoom level is maximum and the tiles are not shown. We should choose a sensible zoom value.

Test against several GDAL versions

CircleCI allows using custom Docker images, which would be extremely useful to test against different versions of GDAL:

https://circleci.com/docs/2.0/custom-images/#overview

Also, the Workflows functionality is more similar to what can be found in GitLab CI:

https://circleci.com/docs/2.0/workflows/

Release 0.2

If you think something should go into v0.2.0, say it here.

Merge of rasters that don't share projection fails

it fails becasue of this line https://github.com/satellogic/telluric/blob/master/telluric/georaster.py#L223

the _prepare_other_raster returns None

Plotting FileCollection with Jupyter Notebook

Hello and thanks for your amazing work!

Is it possible to plot FileCollection with Jupyter Notebook? Now I get the following error:

import os
import telluric as tl

os.environ['GDAL_DATA'] = '/usr/share/gdal/2.2'
fp = tl.FileCollection.open('/home/rykov/download/326.geojson')
fp

---------------------------------------------------------------------------
CPLE_OpenFailedError                      Traceback (most recent call last)
rasterio/_base.pyx in rasterio._base._osr_from_crs()

rasterio/_err.pyx in rasterio._err.exc_wrap_int()

CPLE_OpenFailedError: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

During handling of the above exception, another exception occurred:

CRSError                                  Traceback (most recent call last)
~/sandbox/telluric/env/lib/python3.5/site-packages/IPython/core/formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/plotting.py in _repr_html_(self)
    140         warnings.warn(
    141             "Plotting a limited representation of the data, use the .plot() method for further customization")
--> 142         return simple_plot(self)._repr_html_()
    143 
    144     def plot(self, mp=None, **plot_kwargs):

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/plotting.py in simple_plot(feature, mp, **map_kwargs)
     37 
     38     else:
---> 39         folium.GeoJson(mapping(feature), name='geojson', overlay=True).add_to(mp)
     40         shape = feature.envelope.get_shape(WGS84_CRS)
     41         mp.fit_bounds([shape.bounds[:1:-1], shape.bounds[1::-1]])

~/sandbox/telluric/env/lib/python3.5/site-packages/shapely/geometry/geo.py in mapping(ob)
     80 def mapping(ob):
     81     """Returns a GeoJSON-like mapping"""
---> 82     return ob.__geo_interface__

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/collections.py in __geo_interface__(self)
     97     @property
     98     def __geo_interface__(self):
---> 99         return self.to_record(WGS84_CRS)
    100 
    101     def to_record(self, crs):

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/collections.py in to_record(self, crs)
    102         return {
    103             'type': 'FeatureCollection',
--> 104             'features': [feature.to_record(crs) for feature in self],
    105         }
    106 

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/collections.py in <listcomp>(.0)
    102         return {
    103             'type': 'FeatureCollection',
--> 104             'features': [feature.to_record(crs) for feature in self],
    105         }
    106 

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/features.py in to_record(self, crs)
     45             'type': 'Feature',
     46             'properties': self._attributes,
---> 47             'geometry': self.geometry.to_record(crs),
     48         }
     49 

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/vectors.py in to_record(self, crs)
    324 
    325     def to_record(self, crs):
--> 326         data = mapping(self.get_shape(crs))
    327         if data['type'] == 'LinearRing':
    328             data['type'] = 'Polygon'

~/sandbox/telluric/env/lib/python3.5/site-packages/telluric/vectors.py in get_shape(self, crs)
    337 
    338         """
--> 339         if crs == self.crs:
    340             return self._shape
    341         else:

rasterio/_crs.pyx in rasterio._crs._CRS.__eq__()

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

CRSError: Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

python interpreter crashes on raster.save

creating a simple raster crashes the python interpreter on save.
problem can be reproduced wit::

from affine import Affine
import numpy as np
import telluric as tl
shape = (1, 20, 20)
arr = np.ones(shape)
aff = Affine.scale(1,-1)
rr = tl.GeoRaster2(image=arr, affine=aff, crs=tl.constants.WEB_MERCATOR_CRS)
rr.save('demo.tif')

Properly represent geometries that cross the antimeridian

The first polygon does not cross the antimeridian, the second does:

This might be an intrinsic limitation of the projection, but anyway we can do better than returning a clearly wrong result without warning.

empty_from_roi is not preserving affine

Checkout this code:

empty = tl.GeoRaster2.empty_from_roi(roi=raster.footprint(), resolution=raster.resolution(), band_names=raster.band_names,dtype=raster.dtype)
raster.footprint().almost_equals(empty.footprint())

returns False

I would expect it to return True

Debugging further discovers the following code in rasterize:

minx, miny, maxx, maxy = footprint.get_shape(footprint.crs).bounds
dest_resolution = raster.resolution()
affine = Affine.translation(minx, maxy) * Affine.scale(dest_resolution, -dest_resolution)

but:

affine.almost_equals(raster.affine)

return False

more details:

raster_origin = Point(2, 3)
affine = Affine.translation(raster_origin.x, raster_origin.y)
crs = {'init': 'epsg:32620'}
raster.shape = (1,256,256)

implement GeoRaster2 crs conversion

Required some way to convert a GeoRaster from one crs to another (i.e. raster would be identical: crs & affine modified, but reprents the same coordinates in the world).
e.g. convert raster from mercator to UTM or to Geographic.

Note that there is a project() method, which might be intended for this but it is currently empty. AFAIU the project() method of the old GeoRaster is a generalized version of what I propose here.

First release TODO list

Complete docs: API Reference, (short) User Guide, Development Docs
Add a License ❗
Notebook cleanup
README: Better description, badges
Preparation for pip 10.0 and PEP 518 (pyproject.toml)
Tag v0.1.0rc1, test in https://test.pypi.org/, push v0.1.0 to https://pypi.org
Announcement on python-announce

Left out: wheel auto deployment.

Reproject of rasters loads the raster to memory

Today reprojecting a raster forces to load the raster to memory, it is bad if you want to project a raster that is big

The project code I am referring is:
https://github.com/satellogic/telluric/blob/master/telluric/georaster.py#L945

it is loading everything to memory because of this line:
https://github.com/satellogic/telluric/blob/master/telluric/georaster.py#L969

and this line:
https://github.com/satellogic/telluric/blob/master/telluric/georaster.py#L988

We need to find a way to reproject without loading to memory

A hint I found:
https://rasterio.readthedocs.io/en/latest/topics/reproject.html
we should send a reastio band object instead of the image when image is not loaded

Simplify GeoRaster reproject api

def reproject(self, dest_crs, preserve_pixel_size=True, dtype=None, resampling=Resampling.cubic):
   # If preserve_pixel_size=False, then the resolution is preserved
   if preserve_pixel_size:
       dest_affine = ...
   else:
       dest_affine = ...

   new_width = ...
   new_height = ...

   return self._reproject(new_width, new_height, dest_affine, dtype, dest_crs, resampling)

Retain internal mask while reprojecting

Internal mask is lost after reprojection rasters which are not loaded into memory. Possible ways to solve this solution: https://gist.github.com/drnextgis/ec24f73a895d8602a2d55d2b8d56bd3c

FeatureCollection.save does not retain types

Expected behavior and actual behavior

When saving a FeatureCollection.save and then FileCollection.open, the user would expect to recover exactly the same collection. However, all the resulting types are converted to strings.

Steps to reproduce the problem

In [22]: list(fc)
Out[22]: 
[GeoFeature(Point, {'attr1': 1, 'attr2': None}),
 GeoFeature(Point, {'attr2': '2', 'attr1': None})]

In [23]: fc.save("/tmp/test.shp")

In [24]: fc_saved = tl.FileCollection.open("/tmp/test.shp")

In [25]: list(fc_saved)
Out[25]: 
[GeoFeature(Point, {'attr1': '1', 'attr2': ''}),
 GeoFeature(Point, {'attr1': '', 'attr2': '2'})]

Operating system and installation details

Linux Mint 18.3, latest telluric.

test_geographic_crop fails with rasterio 1.0b1

self = <test_georaster_tiling.GeoRasterCropTest testMethod=test_geographic_crop>

    def test_geographic_crop(self):
        raster = self.geographic_raster()
        rhombus_on_image = Polygon([[0, 2], [1, 1], [2, 2], [1, 3]])  # in pixels
        rhombus_world = raster.to_world(rhombus_on_image)
        cropped = raster.crop(rhombus_world)
        r = raster[0:2, 1:3]
>       assert cropped == r
E       AssertionError: assert <telluric.geo...x7f640559b5c0> == <telluric.geor...x7f640559b630>
E         Full diff:
E         - <telluric.georaster.GeoRaster2 object at 0x7f640559b5c0>
E         ?                                                     ^^
E         + <telluric.georaster.GeoRaster2 object at 0x7f640559b630>
E         ?                                                     ^^

tests/test_georaster_tiling.py:420: AssertionErro

Add support to compute geodetic distances and lengths

Merge to produce metadata of the merge process

We need to know know for each pixel that was produced using Merge, MergeAll for which GeoFeature it came.

The use-case here is to apply to a pixel generated by merge the attributes defined in the GeoFeature the pixel came from

I was thinking maybe to return a raster with integers as values, each integer is the index of the GeoFeature in the FeatureCollection

But we can discuss other ideas as well

Improve plotting

At the moment, default plotting is using a simple Folium map which has bad performance, and more complex .plot using ipyleaflet is still limiting the maximum number of features:

https://github.com/satellogic/telluric/blob/master/telluric/collections.py#L179-L183

We have to re-evaluate if this is needed.

Also, mapboxgl-jupyter just released a new version with Choropleth support:

https://github.com/mapbox/mapboxgl-jupyter/releases/tag/0.6.0

It still does not support all geometry types (mapbox/mapboxgl-jupyter#4) but it's in a better shape than two months ago.

GeoRaster2 does not warn about different mask per band

Expected behavior and actual behavior

rasterio supports mask per dataset and not per band. However, no warning or error is emitted when attempting to create a raster with a different mask per band, which will not retain it on write & open.

Steps to reproduce the problem

In [1]: import telluric as tl
   ...: from telluric.constants import WGS84_CRS
   ...: 
   ...: 

In [2]: import numpy as np

In [3]: from affine import Affine

In [4]: affine = Affine.translation(0, 2) * Affine.scale(1, -1)

In [5]: raster = tl.GeoRaster2(
   ...:     image=np.array([
   ...:         [
   ...:             [100, 200],
   ...:             [100, 200]
   ...:         ],
   ...:         [
   ...:             [110, 0],
   ...:             [110, 0]
   ...:         ]
   ...:     ], dtype=np.uint8),
   ...:     affine=affine,
   ...:     crs=WGS84_CRS,
   ...:     band_names=['red', 'green'],
   ...:     nodata=0
   ...: )

In [6]: print(raster.image)
[[[100 200]
  [100 200]]

 [[110 --]
  [110 --]]]

In [7]: raster.save("/tmp/test_mask.tiff")

In [9]: print(tl.GeoRaster2.open("/tmp/test_mask.tiff").image)
[[[100 --]
  [100 --]]

 [[110 --]
  [110 --]]]

The same happens with rasters created from masked arrays:

In [11]: raster = tl.GeoRaster2(
    ...:     image=np.ma.masked_array([
    ...:         [
    ...:             [100, 200],
    ...:             [100, 200]
    ...:         ],
    ...:         [
    ...:             [110, 0],
    ...:             [110, 0]
    ...:         ]
    ...:     ], [
    ...:         [
    ...:             [False, False],
    ...:             [False, False]
    ...:         ],
    ...:         [
    ...:             [False, True],
    ...:             [False, True]
    ...:         ]
    ...:     ], dtype=np.uint8),
    ...:     affine=affine,
    ...:     crs=WGS84_CRS,
    ...:     band_names=['red', 'green'],
    ...: )

In [12]: print(raster.image)
[[[100 200]
  [100 200]]

 [[110 --]
  [110 --]]]

In [13]: raster.save("/tmp/test_mask2.tiff")

In [14]: print(tl.GeoRaster2.open("/tmp/test_mask2.tiff").image)
[[[100 --]
  [100 --]]

 [[110 --]
  [110 --]]]

Operating system and installation details

Latest telluric, Linux.

Implement proper geodesic buffering

Expected behavior and actual behavior

When doing GeoVector.buffer the user would expect for proper geodesic buffering to happen, as described in this article:

https://www.esri.com/news/arcuser/0111/geodesic.html

Instead, the wrong result is displayed.

Steps to reproduce the problem

tl.GeoVector(
    Point(125, 39.36827914916014), WGS84_CRS
).reproject(WEB_MERCATOR_CRS).buffer(10000e3)

Operating system and installation details

Linux, latest telluric.

crs is sometimes dict, sometimes rasterio.crs.CRS

We are doing internal conversions in some places, but stronger checks are needed, especially now that comparison between CRS and dict has been deprecated in rasterio:

https://github.com/mapbox/rasterio/blob/832b697948c715b9cc04aa329368cb65ad02082e/CHANGES.txt#L4-L30

Reprojecting south-up images

All works as expected for north-up images. But for south-up (does telluric intent to support such images?) we have some problems are appeared at least while reprojecting. I've tried to rescale images and faced with the following issues:

in-memory reprojection switches sign of affine e
result image of non-in-memory reprojection is matched with rasterio warp but not with gdalwarp
result image of in-memory reprojection is not matched neither with rasterio warp nor gdalwarp

I have prepared two notebooks:

north-up image (all works without any issues): https://gist.github.com/drnextgis/e6deadf0f83f30f9d2ea438b3c42f556
south-up image (has issues described above): https://gist.github.com/drnextgis/b3476cb0a102a5be3888e89f9fc5ddab

Such south-up images are used in our tests, for example:

telluric/tests/test_georaster.py

Line 31 in 4090823

    
           some_raster = GeoRaster2(some_image_2d, affine=some_affine, crs=some_crs, band_names=['r'])

Upload Docker image to dockerhub

Getting histogram fails

In [1]: import telluric as tl

In [2]: rr  = tl.GeoRaster2.open('tests/data/raster/rgb.tif')

In [3]: rr.histogram()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-3c91277270c7> in <module>()
----> 1 rr.histogram()

~/sandbox/telluric/telluric/telluric/georaster.py in histogram(self)
   1193 
   1194     def histogram(self):
-> 1195         histogram = {band: self._histogram(band) for band in self.band_names}
   1196         return Histogram(histogram)
   1197 

~/sandbox/telluric/telluric/telluric/georaster.py in <dictcomp>(.0)
   1193 
   1194     def histogram(self):
-> 1195         histogram = {band: self._histogram(band) for band in self.band_names}
   1196         return Histogram(histogram)
   1197 

~/sandbox/telluric/telluric/telluric/georaster.py in _histogram(self, band)
   1204             raise GeoRaster2NotImplementedError('cant calculate histogram for type %s' % self.image.dtype)
   1205 
-> 1206         band_image = self.limit_to_bands(band).image
   1207         return np.histogram(band_image[~band_image.mask], range=(0, length), bins=length)[0]
   1208 

~/sandbox/telluric/telluric/telluric/georaster.py in limit_to_bands(self, bands)
   1082 
   1083     def limit_to_bands(self, bands):
-> 1084         subimage = self.subimage(bands)
   1085         return self.copy_with(image=subimage, band_names=bands)
   1086 

~/sandbox/telluric/telluric/telluric/georaster.py in subimage(self, bands)
    387             bands = bands.split(",")
    388 
--> 389         missing_bands = set(bands) - set(self.band_names)
    390         if missing_bands:
    391             raise GeoRaster2Error('requested bands %s that are not found in raster' % missing_bands)

TypeError: 'int' object is not iterable

Irrelevant exception while saving raster without affine

In [1]: import telluric as tl

In [2]: import numpy as np

In [3]: shape = (1, 20, 20)

In [4]: arr = np.ones(shape)

In [5]: rr = tl.GeoRaster2(image=arr, crs=tl.constants.WEB_MERCATOR_CRS)

In [6]: rr.save('demo.tif')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-caf1c532faea> in <module>()
----> 1 rr.save('demo.tif')

~/sandbox/telluric/telluric/telluric/georaster.py in save(self, filename, tags, **kwargs)
    570 
    571                 params = {
--> 572                     'mode': "w", 'transform': self.affine, 'crs': self.crs,
    573                     'driver': driver, 'width': size[2], 'height': size[1], 'count': size[0],
    574                     'dtype': dtype_map[self.image.dtype.type],

~/sandbox/telluric/telluric/telluric/georaster.py in affine(self)
    488         """Raster affine."""
    489         if self._affine is None:
--> 490             self._populate_from_rasterio_object(read_image=False)
    491         return self._affine
    492 

~/sandbox/telluric/telluric/telluric/georaster.py in _populate_from_rasterio_object(self, read_image)
    434 
    435     def _populate_from_rasterio_object(self, read_image):
--> 436         with self._raster_opener(self._filename) as raster:  # type: rasterio.DatasetReader
    437             self._affine = copy(raster.transform)
    438             self._crs = copy(raster.crs)

~/sandbox/telluric/telluric/telluric/georaster.py in _raster_opener(cls, filename, *args, **kwargs)
    406             warnings.simplefilter('ignore', UserWarning)
    407             try:
--> 408                 return rasterio.open(filename, *args, **kwargs)
    409             except (rasterio.errors.RasterioIOError, rasterio._err.CPLE_BaseError) as e:
    410                 raise GeoRaster2IOError(e)

~/sandbox/telluric/env/lib/python3.5/site-packages/rasterio/__init__.py in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, **kwargs)
    156     if not isinstance(fp, string_types):
    157         if not (hasattr(fp, 'read') or hasattr(fp, 'write')):
--> 158             raise TypeError("invalid path or file: {0!r}".format(fp))
    159     if mode and not isinstance(mode, string_types):
    160         raise TypeError("invalid mode: {0!r}".format(mode))

TypeError: invalid path or file: None

Non shrunk masks are assumed in some parts of the code

There is one place where we take care of a shrunk mask (that is: when numpy.ma decides to replace an array of False values with a np.ma.nomask scalar value):

telluric/telluric/georaster.py

Lines 1096 to 1097 in 7b393a0

    
           if mask is np.ma.nomask: 
        
               mask = np.zeros_like(self.image.data, dtype=bool)

However, there are other places where we don't protect ourselves against this case, for example:

telluric/telluric/georaster.py

Line 218 in 7b393a0

other_values_mask = (other_image.mask[0] | (~one.image.mask[0]))

We should treat this issue in a more consistent way. We have already noticed one case where a merge_all fails with properly formatted data.

an empty feature collection fails on `property_names`

whenever you try to access property_name() of an empty FeatureCollection you get an exception: KeyError: 'pop from an empty set'

reproduce:

import telluric as tl
fc = tl.FeatureCollection([])
fc.property_names  # this fails with KeyError

Loading and saving of a GeoRaster does not preserve raster attributes

There are cases where attributes of a loaded raster are not preserve on saving:

When loading a raster without tiling and /or overviews, the raster is saved tiled with overviews.
When the overviews of the original raster are different than the on the save method

The source of that is we don't record this raster attributes on loading

Affine is not saved on when affine equals (1.00, 0.00, 0.00, 0.00,-1.00, 0.00)

When saving a GeoRaster2 with affine:

| 1.00, 0.00, 0.00|
| 0.00,-1.00, 0.00|
| 0.00, 0.00, 1.00|

reading the saved raster returns affine:

| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

Sample code:

from affine import Affine
import numpy as np
import telluric as tl
shape = (1, 500, 500)
arr = np.ones(shape)
aff = Affine.scale(1,-1)
raster1 = tl.GeoRaster2(image=arr, affine=aff, crs=tl.constants.WEB_MERCATOR_CRS)
raster1.save('demo.tif')

raster2 = tl.GeoRaster2.open('demo.tif')

raster1.affine == raster2.affine

Use maximum width or height for calculating overviews factors

It is common practice to calculate overviews factors using maximum width or heigh of raster. But in _overviews_factors method here min function is used. Was it done intentionally?

`get_tiles` occasionally returns white borders around the tile

see attached example of tiles generated using get_tile and presented in Folium:

Improve GeoRaster dtype conversion and stretching

At the moment, conversions to/from float32/64 are not implemented. This is partly because the maximum values of just single precision floats are too high (~10^38), and therefore doing a proper stretching is difficult.

On the other hand, we usually face the problem of wanting to stretch a sequence of images with the same values range.

Therefore, it should be possible to specify a range when stretching/changing the dtype of a GeoRaster2.

Merging non overlapping rasters is order dependant

Expected behavior and actual behavior

I expected that merging non overlapping rasters would always give the same result regardless of order, but it's not the case.

Steps to reproduce the problem

affine = Affine.translation(0, 2) * Affine.scale(1, -1)

rs1 = tl.GeoRaster2(
    image=np.array([[
        [100, 0],
        [100, 0]
    ]], dtype=np.uint8),
    affine=affine,
    crs=WGS84_CRS,
    band_names=['red'],
)

rs2 = tl.GeoRaster2(
    image=np.array([[
        [110, 0],
        [110, 0]
    ]], dtype=np.uint8),
    affine=affine,
    crs=WGS84_CRS,
    band_names=['green'],
)

rs3 = tl.GeoRaster2(
    image=np.array([[
        [0, 200],
        [0, 200]
    ]], dtype=np.uint8),
    affine=affine,
    crs=WGS84_CRS,
    band_names=['red'],
)

rs4 = tl.GeoRaster2(
    image=np.array([[
        [0, 210],
        [0, 210]
    ]], dtype=np.uint8),
    affine=affine,
    crs=WGS84_CRS,
    band_names=['green'],
)

print(tl.georaster.merge_all([rs1, rs2, rs3, rs4], roi).image)
print(tl.georaster.merge_all([rs1, rs3, rs2, rs4], roi).image)

Operating system and installation details

Latest telluric, Linux.

FeatureCollection.save has side effects on heterogeneous collections

When a FeatureCollection contains features with different attributes, .save "fills" the missing ones with None:

In [17]: import telluric as tl

In [18]: from shapely.geometry import Point

In [19]: fc = tl.FeatureCollection([
    ...: tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr1': 1}),
    ...: tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr2': '2'})
    ...: ])

In [20]: list(fc)
Out[20]: [GeoFeature(Point, {'attr1': 1}), GeoFeature(Point, {'attr2': '2'})]

In [21]: fc.save("/tmp/test.shp")

In [22]: list(fc)
Out[22]: 
[GeoFeature(Point, {'attr1': 1, 'attr2': None}),
 GeoFeature(Point, {'attr2': '2', 'attr1': None})]

transform_attributes function doesn't respect None values

Steps to reproduce the problem

>>> import telluric as tl
>>> from shapely.geometry import Point
>>> fc = tl.FeatureCollection([tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr1': '2018-05-19'}), tl.GeoFeature(tl.GeoVector(Point(0, 0)), {'attr1': None})])
>>> fc.save('/tmp/test3.geojson')
Warning 6: dataset /tmp/test3.geojson does not support layer creation option ENCODING
>>> fc_saved = tl.FileCollection.open("/tmp/test3.geojson")
>>> list(fc_saved)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/rykov/sandbox/telluric/telluric/telluric/collections.py", line 418, in __iter__
    yield GeoFeature.from_record(record, source.crs, source.schema)
  File "/home/rykov/sandbox/telluric/telluric/telluric/features.py", line 101, in from_record
    attributes = transform_attributes(record["properties"], schema)
  File "/home/rykov/sandbox/telluric/telluric/telluric/features.py", line 33, in transform_attributes
    new_attributes[attr_name] = parse_date(attr_value).date()
  File "/home/rykov/sandbox/telluric/env/lib/python3.5/site-packages/dateutil/parser.py", line 1182, in parse
    return DEFAULTPARSER.parse(timestr, **kwargs)
  File "/home/rykov/sandbox/telluric/env/lib/python3.5/site-packages/dateutil/parser.py", line 556, in parse
    res, skipped_tokens = self._parse(timestr, **kwargs)
  File "/home/rykov/sandbox/telluric/env/lib/python3.5/site-packages/dateutil/parser.py", line 675, in _parse
    l = _timelex.split(timestr)         # Splits the timestr into tokens
  File "/home/rykov/sandbox/telluric/env/lib/python3.5/site-packages/dateutil/parser.py", line 192, in split
    return list(cls(s))
  File "/home/rykov/sandbox/telluric/env/lib/python3.5/site-packages/dateutil/parser.py", line 61, in __init__
    '{itype}'.format(itype=instream.__class__.__name__))
TypeError: Parser must be a string or character stream, not NoneTyp

Operating system and installation details

Linux, latest telluric

Cropping raster with GeoVector bigger than world bounds fails

Expected behavior and actual behavior

Cropping any raster with any GeoVector should work. Instead, if the GeoVector is bigger than the world, an error is raised.

simple_world.geojson.zip

Steps to reproduce the problem

In [1]: import telluric as tl

In [2]: rs_simple = tl.GeoRaster2.open("./tests/data/raster/creaf_gmap.tif")

In [3]: gv = tl.GeoVector.from_geojson("/tmp/simple_world.geojson")

In [4]: rs_simple.crop(gv)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-53bf33faf479> in <module>()
----> 1 rs_simple.crop(gv)

~/Development/telluric/telluric/georaster.py in crop(self, vector, resolution)
    710         :return: GeoRaster
    711         """
--> 712         bounds, window = self._vector_to_raster_bounds(vector)
    713         if resolution:
    714             xsize, ysize = self._resolution_to_output_shape(bounds, resolution)

~/Development/telluric/telluric/georaster.py in _vector_to_raster_bounds(self, vector, boundless)
    720     def _vector_to_raster_bounds(self, vector, boundless=False):
    721         # bounds = tuple(round(bb) for bb in self.to_raster(vector).bounds)
--> 722         window = self.window(*vector.get_shape(self.crs).bounds).round_offsets().round_shape(op='ceil')
    723         (ymin, ymax), (xmin, xmax) = window.toranges()
    724         bounds = (xmin, ymin, xmax, ymax)

~/.miniconda36/envs/telluric36/lib/python3.6/site-packages/rasterio/windows.py in round_offsets(self, op, pixel_precision)
    708             raise WindowError("operator must be 'ceil' or 'floor'")
    709         else:
--> 710             return Window(operator(round(self.col_off, pixel_precision)),
    711                           operator(round(self.row_off, pixel_precision)),
    712                           self.width, self.height)

ValueError: cannot convert float NaN to integer

Notice that with a simple .buffer(-1) the problem is solved in this case:

In [5]: gv.envelope
Out[5]: GeoVector(shape=POLYGON ((-180 -90, 180 -90, 180 83.70384706880299, -180 83.70384706880299, -180 -90)), crs=CRS({'init': 'epsg:4326'}))

In [6]: gv.buffer(-1).envelope  # slow
Out[6]: GeoVector(shape=POLYGON ((-179 -89, 179 -89, 179 82.70053994948313, -179 82.70053994948313, -179 -89)), crs=CRS({'init': 'epsg:4326'}))

In [7]: rs_simple.crop(gv.buffer(-1))
Out[7]: <telluric.georaster.GeoRaster2 at 0x7eff872bf630>

Operating system and installation details

Linux Mint 18.3, latest telluric.

Crop size incorrectly computed for WGS84

Expected behavior and actual behavior

I expected raster.crop(roi) to produce the same result as raster.crop(roi, raster.resolution()), but when raster.crs is WGS84_CRS, the size is incorrectly computed, which might result in a MemoryError.

Steps to reproduce the problem

In [1]: import numpy as np
   ...: 
   ...: import telluric as tl
   ...: from telluric.constants import WGS84_CRS
   ...: 
   ...: 

In [2]: roi = tl.GeoVector.from_bounds(xmin=0, ymin=0, xmax=2, ymax=2, crs=WGS84_CRS)
   ...: resolution = 1.0
   ...: 
   ...: 

In [4]: rs0 = tl.GeoRaster2.empty_from_roi(roi, resolution)

In [6]: rs0.crop(roi).image
Out[6]: 
masked_array(
  data=[[[--, --],
         [--, --]]],
  mask=[[[ True,  True],
         [ True,  True]]],
  fill_value=999999,
  dtype=uint8)

In [7]: rs0.crop(roi, resolution).image
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-7-0b900c87a084> in <module>()
----> 1 rs0.crop(roi, resolution).image

~/Development/telluric/telluric/georaster.py in crop(self, vector, resolution)
    773             xsize, ysize = (None, None)
    774 
--> 775         return self.pixel_crop(bounds, xsize, ysize, window=window)
    776 
    777     def _vector_to_raster_bounds(self, vector, boundless=False):

~/Development/telluric/telluric/georaster.py in pixel_crop(self, bounds, xsize, ysize, window)
    828         """
    829         if self._image is not None:
--> 830             return self._crop(bounds, xsize=xsize, ysize=ysize)
    831         else:
    832             window = window or rasterio.windows.Window(bounds[0],

~/Development/telluric/telluric/georaster.py in _crop(self, bounds, xsize, ysize)
    851         if xsize and ysize:
    852             if not (xsize == out_raster.width and ysize == out_raster.height):
--> 853                 out_raster = out_raster.resize(dest_width=xsize, dest_height=ysize)
    854         return out_raster
    855 

~/Development/telluric/telluric/georaster.py in resize(self, ratio, ratio_x, ratio_y, dest_width, dest_height, dest_resolution, resampling)
    908             ratio_x, ratio_y = ratio, ratio
    909 
--> 910         return self._resize(ratio_x, ratio_y, resampling)
    911 
    912     def _resize(self, ratio_x, ratio_y, resampling):

~/Development/telluric/telluric/georaster.py in _resize(self, ratio_x, ratio_y, resampling)
    915         new_height = int(np.ceil(self.height * ratio_y))
    916         dest_affine = self.affine * Affine.scale(1 / ratio_x, 1 / ratio_y)
--> 917         return self.reproject(new_width, new_height, dest_affine, resampling=resampling)
    918 
    919     def to_pillow_image(self, return_mask=False):

~/Development/telluric/telluric/georaster.py in reproject(self, new_width, new_height, dest_affine, dtype, dst_crs, resampling)
    951         dtype = dtype or self.image.data.dtype
    952         dest_image = np.ma.masked_array(
--> 953             data=np.empty([self.num_bands, new_height, new_width], dtype=np.float32),
    954             mask=np.empty([self.num_bands, new_height, new_width], dtype=bool)
    955         )

MemoryError:

The reason is this:

In [11]: bounds, window = rs0._vector_to_raster_bounds(roi)

In [12]: bounds, window
Out[12]: ((0, 0, 2, 2), Window(col_off=0, row_off=0, width=2, height=2))

In [13]: rs0._resolution_to_output_shape(bounds, resolution)
Out[13]: (222222.0, 222222)

Operating system and installation details

Linux, latest telluric

Meaning of _calc_overviews_factors's blocksize argument

@arielze I have a question related to the algorithm of overviews factors calculation:

telluric/telluric/util/raster_utils.py

Line 11 in d7206c3

def _calc_overviews_factors(one, blocksize=256):

What is a blocksize argument? Is it maximum width or height of the smallest overview level (as done here, see minsize option). If it is than we should modify _calc_overviews_factors a little bit, because it returns unexpected results. For example:

In [1]: import numpy as np

In [2]: import telluric as tl

In [3]: shape = (1, 1024, 256)

In [4]: arr = np.ones(shape)

In [5]: img = tl.GeoRaster2(image=arr, crs=tl.constants.WEB_MERCATOR_CRS)

In [6]: img._overviews_factors(blocksize=256)
Out[6]: [2, 4, 8]

I've expected to get [2, 4] here. Example of gdaladdo results:

$ gdalinfo --version
GDAL 2.3.0, released 2018/05/04

$ gdaladdo demo.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo demo.tif
Driver: GTiff/GeoTIFF
Files: demo.tif
Size is 256, 1024
Coordinate System is:
...
Metadata:
  AREA_OR_POINT=Area
  telluric_band_names=[0]
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1024.0)
Upper Right (  256.0,    0.0)
Lower Right (  256.0, 1024.0)
Center      (  128.0,  512.0)
Band 1 Block=256x256 Type=Float64, ColorInterp=Gray
  Overviews: 128x512, 64x256
  Mask Flags: PER_DATASET 
  Overviews of mask band: 128x512, 64x256

	if mask is np.ma.nomask:
	mask = np.zeros_like(self.image.data, dtype=bool)

satellogic / telluric Goto Github PK

telluric's Introduction

telluric

Overview

Installation

Development

Support

telluric's People

Contributors

Stargazers

Watchers

Forkers

telluric's Issues

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Expected behavior and actual behavior

Steps to reproduce the problem

Operating system and installation details

Recommend Projects

Recommend Topics

Recommend Org