search

OpenGeoscience / geonotebook

A Jupyter notebook extension for geospatial visualization and analysis
Apache License 2.0
1.08k stars 141 forks source link

Getting raster data from timeseries rasters under annotation layer #128

Open digitaltopo opened 7 years ago

digitaltopo commented 7 years ago

I have a layer created using a raster data collection RasterDataCollection() that I'm using to create a time series, following the 03_Time_Series.ipynb example notebook.

I'm trying to grab raster data that is bounded by an annotation for each raster included in the time series (each date in timeseries). In the 02_Raster_Data.ipynb example notebook, there is an example of getting raster data from the annotation using:

layer, data = next(M.layers.annotation.rectangles[0].data)
   data

I'm seeing this output:

masked_array(data =
 [],
             mask =
 False,
       fill_value = -3.40282e+38)

But not sure where to go from here. How does the annotation layer relate to data in other layers? Thanks!

kotfic commented 7 years ago

@digitaltopo Assuming you've set up your RasterDataCollection correctly then next(M.layers.annotation.rectangles[0].data) should return data from the first rectangle annotation from the first layer. If you've only got one layer on the map, and that layer was created from the RasterDataCollection, and your rectangle annotation is actually over an area that has meaningful data, then the first dimension of the masked array should be each of the time series slices sub-selected by the annotation.

Would it be possible to provide a minimum working example of the error you're seeing?

digitaltopo commented 7 years ago

Thanks, this was my assumption too, but the masked_array seems to be empty. Here is the basic code I'm using to create:

# Add layer with data from rasters
DATA_DIR="data/"
tifs = sorted(glob(DATA_DIR + '/**/*.tif'))
rdc = RasterDataCollection(tifs)
M.add_layer(rdc[:,[1]], 'TimeSeriesLayer', colormap=cmap, opacity=0.8, gamma=1.0)

# Time series control
def render_timeseries(idx=0):
  M.layers[0].idx(idx)

interact(render_timeseries, idx=(0, len(M.layers[0].data) - 1))

layer, data = next(M.layers.annotation.rectangles[0].data)
layer

The "error" is that the masked_array seems to not include any data. 

masked_array(data =
 [],
             mask =
 False,
       fill_value = -3.40282e+38)

Is this the correct process, or perhaps there is an issue with my raster? The rasters display fine on the map etc.

kotfic commented 7 years ago

@digitaltopo No, that looks right to me. If it's not too much trouble could you try going into tests/integration folder and running ./get_test_data.sh? This will download a couple of geotiffs into the data/ directory.

Then you should be able to run jupyter notebook and select the TimeSeries.ipynb test notebook. Run the first two cells (and make sure your output looks like the image that is displayed in the test).

Then:

  1. make a new cell
  2. select a region of the data with the rectangle annotation tool
  3. run next(M.layers.annotation.rectangles[0].data)
  4. Check that the output contains a masked_array with data

Something like:

(<TimeSeriesLayer('WELD')>, masked_array(data =
  [[[ 0.0872      0.0384      0.0507     ...,  0.14309999  0.0607      0.0384    ]
   [ 0.2139      0.2674      0.1585     ...,  0.1089      0.0557      0.1585    ]
   [ 0.25749999  0.20810001  0.1534     ..., -0.0218      0.0319      0.0512    ]
   ..., 
   [ 0.3725      0.2403      0.0798     ...,  0.3646      0.7403      0.6559    ]
   [ 0.1189      0.0849      0.0491     ...,  0.6397      0.55930001
     0.47600001]
   [ 0.0864      0.15270001  0.2102     ...,  0.66549999  0.66469997
     0.56199998]]

  [[ 0.2581      0.16590001  0.088      ...,  0.1112      0.0831      0.0165    ]
   [ 0.0883      0.2938      0.088      ...,  0.0452      0.0398      0.14489999]
   [ 0.16769999  0.1673      0.1673     ...,  0.0552     -0.0009      0.1278    ]
   ..., 
   [ 0.57999998  0.31        0.28389999 ...,  0.30939999  0.61220002
     0.60030001]
   [ 0.22939999  0.12809999  0.0864     ...,  0.63150001  0.55510002  0.5043    ]
   [ 0.18970001  0.0643      0.1426     ...,  0.486       0.62959999
     0.52899998]]

  [[ 0.4382      0.0772      0.1224     ...,  0.0762      0.0894      0.1031    ]
   [ 0.1313      0.1013      0.0524     ...,  0.0883      0.1033      0.13330001]
   [ 0.0829      0.1224      0.1004     ...,  0.0287      0.0341      0.1114    ]
   ..., 
   [ 0.49439999  0.31619999  0.3488     ...,  0.1875      0.382       0.42809999]
   [ 0.26949999  0.0615      0.1996     ...,  0.46900001  0.26710001
     0.34209999]
   [ 0.1195      0.0165      0.28259999 ...,  0.2217      0.49090001
     0.25830001]]],
              mask =
  False,
        fill_value = -9999.0))

If this works then it might be some incompatibility with the files your using in which case i'd love to see if I can reproduce the error locally.

digitaltopo commented 7 years ago

Yes it seems to work as you've posted with the test rasters. What property on the TimeSeriesLayer is it looking at, maybe I can try to see how my rasters differ from the test rasters. Thanks!

dorukozturk commented 6 years ago

@digitaltopo We merged a PR which should solve this issue. Can you pull the latest geonotebook master and try with your data now?

digitaltopo commented 6 years ago

@dorukozturk I pulled master and reinstalled but I am now getting an error when I start geonotebook: 

ImportError: /home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/rasterio/_warp.so: undefined symbol: _ZTINSt8ios_base7failureB5cxx11E

Full stack trace:

[W 10:22:38.135 NotebookApp] server_extensions is deprecated, use nbserver_extensions
[W 10:22:38.811 NotebookApp] Error loading server extension geonotebook
    Traceback (most recent call last):
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/notebook/notebookapp.py", line 1271, in init_server_extensions
        mod = importlib.import_module(modulename)
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/importlib/__init__.py", line 37, in import_module
        __import__(name)
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/__init__.py", line 6, in <module>
        from .config import Config
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/config.py", line 79, in <module>
        Config.register_vis_server(ep.name, ep.load())
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/pkg_resources/__init__.py", line 2306, in load
        return self.resolve()
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/pkg_resources/__init__.py", line 2312, in resolve
        module = __import__(self.module_name, fromlist=['__name__'], level=0)
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/vis/__init__.py", line 1, in <module>
        from .geoserver import Geoserver
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/vis/geoserver/__init__.py", line 1, in <module>
        from .geoserver import Geoserver
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/vis/geoserver/geoserver.py", line 4, in <module>
        from geonotebook.wrappers import RasterData
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/wrappers/__init__.py", line 1, in <module>
        from .raster import RasterData, RasterDataCollection
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/wrappers/raster.py", line 177, in <module>
        RasterData.discover_concrete_types()
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/wrappers/raster.py", line 29, in discover_concrete_types
        cls.register(ep.name, ep.load())
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/pkg_resources/__init__.py", line 2306, in load
        return self.resolve()
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/pkg_resources/__init__.py", line 2312, in resolve
        module = __import__(self.module_name, fromlist=['__name__'], level=0)
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/wrappers/file_reader.py", line 11, in <module>
        from geonotebook.utils import transform_coordinates
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/geonotebook/utils.py", line 4, in <module>
        from rasterio.warp import transform
      File "/home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/rasterio/warp.py", line 13, in <module>
        from rasterio._warp import (
    ImportError: /home/USERDIR/.applications/anaconda3/envs/geo2/lib/python2.7/site-packages/rasterio/_warp.so: undefined symbol: _ZTINSt8ios_base7failureB5cxx11E

My update process was to pull master, run pip install . --upgrade --force-reinstall and run npm install && npm build for the js dir. I also reinstalled rasterio but same issue.

dorukozturk commented 6 years ago

What is your gdal version? You can do:

gdal-config --version

for gdal system version and 

pip freeze

for your gdal-python binding version. The result of gdal-config --version and gdal version on pip freeze should match and should be >= 1.9.

digitaltopo commented 6 years ago

gdal-config --version : 2.1.0 pip freeze: GDAL==2.1.0

dorukozturk commented 6 years ago

I see you are on Conda. Can you try to install rasterio with:

conda install -c conda-forge rasterio