satellogic / telluric Goto Github PK
View Code? Open in Web Editor NEWtelluric is a Python library to manage vector and raster geospatial data in an interactive and easy way
License: MIT License
telluric is a Python library to manage vector and raster geospatial data in an interactive and easy way
License: MIT License
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
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
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.
tl.GeoVector(
Point(125, 39.36827914916014), WGS84_CRS
).reproject(WEB_MERCATOR_CRS).buffer(10000e3)
Linux, latest telluric.
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)
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.
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})]
>>> 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
Linux, latest telluric
Slicing a FeatureCollection
should return another FeatureCollection, but now it returns a list.
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})]
Linux, latest telluric.
pyproject.toml
)Left out: wheel auto deployment.
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:
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.
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'})]
Linux Mint 18.3, latest telluric.
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.
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 --]]]
Latest telluric, Linux.
If you think something should go into v0.2.0, say it here.
it fails becasue of this line https://github.com/satellogic/telluric/blob/master/telluric/georaster.py#L223
the _prepare_other_raster returns None
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.
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?
There are cases where attributes of a loaded raster are not preserve on saving:
The source of that is we don't record this raster attributes on loading
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.
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
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
We should have debug mode that will help us understand what is going on.
In this mode we should see:
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:
e
rasterio warp
but not with gdalwarp
rasterio warp
nor gdalwarp
I have prepared two notebooks:
Such south-up images are used in our tests, for example:
telluric/tests/test_georaster.py
Line 31 in 4090823
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)
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
.
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)
Linux, latest telluric
There are cases when two rasters have the same information, but comparing them returns False
, for example:
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
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
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.
I expected that merging non overlapping rasters would always give the same result regardless of order, but it's not the case.
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)
Latest telluric, Linux.
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
.
On the readthedocs page: http://telluric.readthedocs.io/en/v0.1.0/User%20Guide.html
The comment on In [6]: gv1.area # Real area in meters
this should read In [6]: gv1.area # Real area in square meters
Just read the documents.
Any & All.
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
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
@arielze I have a question related to the algorithm of overviews factors calculation:
telluric/telluric/util/raster_utils.py
Line 11 in d7206c3
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
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
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
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')
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
Cropping any raster with any GeoVector
should work. Instead, if the GeoVector
is bigger than the world, an error is raised.
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>
Linux Mint 18.3, latest telluric.
Ideas after telluric 0.2.0 is released (#102):
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
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
However, there are other places where we don't protect ourselves against this case, for example:
telluric/telluric/georaster.py
Line 218 in 7b393a0
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.
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
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.
Documentation: https://conda-forge.org/docs/recipe.html
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.