Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Drawing coastlines/borders with save_datasets #464

Closed
hany opened this issue Oct 16, 2018 · 26 comments · Fixed by #469
Closed

Drawing coastlines/borders with save_datasets #464

hany opened this issue Oct 16, 2018 · 26 comments · Fixed by #469

Comments

@hany
Copy link

hany commented Oct 16, 2018

Hi there, I'm looking for a way to draw coastlines and political borders with save_datasets(). I see references in other issues to a coast_dir, but not sure what should go in this directory. Cannot find much in the way of documentation either. Any help would be appreciated!

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

You have two options:

  1. You can specify options for the pycoast package in the overlay parameter that is passed to save_datasets. This will be a dictionary of parameters. For available options see https://satpy.readthedocs.io/en/latest/api/satpy.writers.html#satpy.writers.add_overlay (although I'm noticing now that this is terribly formatted for some reason). And for what they do and available values see: https://pycoast.readthedocs.io/en/latest/
  2. Use cartopy to make plots using matplotlib: https://github.com/pytroll/pytroll-examples/blob/master/satpy/Cartopy%20Plot.ipynb

It would be great if you could contribute any documentation to add to satpy to improve things for the various issues you've been making.

@hany
Copy link
Author

hany commented Oct 16, 2018

Thanks @djhoese. I found the docs on overlay here (properly formatted):

"""Add coastline and political borders to image, using *color* (tuple
of integers between 0 and 255).
Warning: Loses the masks !
*resolution* is chosen automatically if None (default), otherwise it should be one of:
+-----+-------------------------+---------+
| 'f' | Full resolution | 0.04 km |
| 'h' | High resolution | 0.2 km |
| 'i' | Intermediate resolution | 1.0 km |
| 'l' | Low resolution | 5.0 km |
| 'c' | Crude resolution | 25 km |
+-----+-------------------------+---------+
"""
. However it's not clear what should go in coast_dir. Going to try digging more in the source, but any help would be appreciated.

It would be great if you could contribute any documentation to add to satpy to improve things for the various issues you've been making.

Would be happy to. Once I get my head wrapped around this (and dependent projects), I'll look into contributing docs.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Ah ok, it is equivalent to the GSHHS_DATA_ROOT mentioned in the pycoast docs.

@hany
Copy link
Author

hany commented Oct 16, 2018

Perfect, thank you. I downloaded the GSHHS and WDBII Shapefiles from
http://www.soest.hawaii.edu/pwessel/gshhs/index.html as specified:

new_scn.save_datasets(overlay={'coast_dir': '/tmp/shapefiles/', 'color': (0, 0, 0)})

I got much farther this time, however I ran into the following exception:

Traceback (most recent call last):
  File "testsatpy.py", line 34, in <module>
    new_scn.save_datasets(overlay={'coast_dir': '/tmp/shapefiles/', 'color': (0, 0, 0)})
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/scene.py", line 1071, in save_datasets
    return writer.save_datasets(datasets, compute=compute, **save_kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 462, in save_datasets
    res = self.save_dataset(ds, compute=False, **kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 563, in save_dataset
    decorate=decorate, fill_value=fill_value)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 354, in get_enhanced_image
    add_overlay(img, dataset.attrs['area'], fill_value=fill_value, **overlay)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 230, in add_overlay
    resolution=resolution, width=width, level=level_coast)
TypeError: add_coastlines() got an unexpected keyword argument 'width'

Looking at pycoast's source, the issue appears to be here:

https://github.com/pytroll/pycoast/blob/master/pycoast/cw_pil.py#L292-L293

As width is not an argument. It is specified here though:

https://github.com/pytroll/pycoast/blob/master/pycoast/cw_agg.py#L381-L383

Traced it back to this:

https://github.com/pytroll/pycoast/blob/master/pycoast/__init__.py#L8-L11

Looks like the import is failing due to the lack of the aggdraw package. After installing it, I was able to get farther until this exception:

Traceback (most recent call last):
  File "testsatpy.py", line 34, in <module>
    new_scn.save_datasets(overlay={'coast_dir': '/Users/hany/hurricane/shapefiles/', 'color': (0, 0, 0)})
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/scene.py", line 1071, in save_datasets
    return writer.save_datasets(datasets, compute=compute, **save_kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 462, in save_datasets
    res = self.save_dataset(ds, compute=False, **kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 565, in save_dataset
    fill_value=fill_value, **kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/geotiff.py", line 192, in save_image
    filename = filename or self.get_filename(**img.data.attrs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 427, in get_filename
    return self.filename_parser.compose(kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/trollsift/parser.py", line 55, in compose
    return compose(self.fmt, keyvals)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/trollsift/parser.py", line 355, in compose
    return formatter.format(fmt, **keyvals)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 190, in format
    return self.vformat(format_string, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 194, in vformat
    result, _ = self._vformat(format_string, args, kwargs, used_args, 2)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 234, in _vformat
    obj, arg_used = self.get_field(field_name, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 299, in get_field
    obj = self.get_value(first, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 256, in get_value
    return kwargs[key]
KeyError: 'name'

Seemed like an error in string formatting for the resultant filename. Manually specifying the filename got me to the end:

new_scn.save_datasets(filename='test.tif', overlay={'coast_dir': '/tmp/shapefiles/', 'color': (0, 0, 0)})
INFO:satpy.writers:Add coastlines and political borders to image.
DEBUG:satpy.writers:Automagically choose resolution i
...
DEBUG:rasterio._base:Dataset <closed DatasetWriter name='test.tif' mode='w'> has been closed.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

What data were you trying to save? The default filename format uses the attrs for name and start_time.

@hany
Copy link
Author

hany commented Oct 16, 2018

I'm using GOES-16 full disk, channel 10. If I do not specify the overlay argument, it renders just fine:

new_scn.save_datasets()
DEBUG:rasterio._base:Dataset <closed DatasetWriter name='C10_20181007_160038.tif' mode='w'> has been closed.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Ah well then that must be a bug. I'll look in to it.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

@hany Are you sure you are doing this with regular old C10? See #449

@hany
Copy link
Author

hany commented Oct 16, 2018

I believe I am. Just downloaded the nc files from the AWS S3 buckets for GOES-16:

DEBUG:satpy.readers.yaml_reader:Assigning to abi_l1b: ['./OR_ABI-L1b-RadF-M3C01_G16_s20182801630380_e20182801641147_c20182801641189.nc', './OR_ABI-L1b-RadF-M3C01_G16_s20182801600380_e20182801611147_c20182801611190.nc', './OR_ABI-L1b-RadF-M3C02_G16_s20182801630380_e20182801641147_c20182801641182.nc', './OR_ABI-L1b-RadF-M3C02_G16_s20182801600380_e20182801611147_c20182801611183.nc', './OR_ABI-L1b-RadF-M3C03_G16_s20182801630380_e20182801641147_c20182801641190.nc', './OR_ABI-L1b-RadF-M3C03_G16_s20182801600380_e20182801611147_c20182801611192.nc', './OR_ABI-L1b-RadF-M3C04_G16_s20182801600380_e20182801611147_c20182801611168.nc', './OR_ABI-L1b-RadF-M3C05_G16_s20182801600380_e20182801611147_c20182801611191.nc', './OR_ABI-L1b-RadF-M3C06_G16_s20182801600380_e20182801611153_c20182801611182.nc', './OR_ABI-L1b-RadF-M3C07_G16_s20182801600380_e20182801611159_c20182801611192.nc', './OR_ABI-L1b-RadF-M3C08_G16_s20182801600380_e20182801611147_c20182801611193.nc', './OR_ABI-L1b-RadF-M3C09_G16_s20182801600380_e20182801611153_c20182801611214.nc', './OR_ABI-L1b-RadF-M3C10_G16_s20182801600380_e20182801611159_c20182801611213.nc', './OR_ABI-L1b-RadF-M3C11_G16_s20182801600380_e20182801611147_c20182801611210.nc', './OR_ABI-L1b-RadF-M3C12_G16_s20182801600380_e20182801611152_c20182801611217.nc', './OR_ABI-L1b-RadF-M3C13_G16_s20182801600380_e20182801611158_c20182801611219.nc', './OR_ABI-L1b-RadF-M3C14_G16_s20182801600380_e20182801611147_c20182801611216.nc', './OR_ABI-L1b-RadF-M3C15_G16_s20182801600380_e20182801611152_c20182801611219.nc', './OR_ABI-L1b-RadF-M3C16_G16_s20182801600380_e20182801611159_c20182801611218.nc']
DEBUG:satpy.composites:Looking for composites config file abi.yaml
DEBUG:satpy.composites:Looking for composites config file visir.yaml
DEBUG:satpy.readers.abi_l1b:Reading in get_dataset C10.
DEBUG:satpy.readers.abi_l1b:Calibrating to brightness temperatures
DEBUG:satpy.scene:Resampling DatasetID(name='C10', wavelength=(7.24, 7.34, 7.44), resolution=2000, polarization=None, calibration='brightness_temperature', level=None, modifiers=())

Is there a better way to verify this?

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Can I see your code from Scene creation to save_datasets?

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

And can you do print(scn['C10']) right before saving?

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Side question: How do you download those files from AWS without an AWS client?

@hany
Copy link
Author

hany commented Oct 16, 2018

Yep, very simple (and messy):

import os
from satpy import Scene
from glob import glob
from dask.diagnostics import ProgressBar

logging.basicConfig()
logger = logging.getLogger()
logger.setLevel(logging.DEBUG)

BASE_DIR = '.'
all_filenames = [glob(fn.replace('C01', 'C[0-9][0-9]*')[:len(BASE_DIR) + 50] + '*.nc') for fn in sorted(glob(os.path.join(BASE_DIR, 'OR*-RadF-*C*.nc')))]
scenes = [Scene(reader='abi_l1b', filenames=filenames) for filenames in all_filenames]
print("Number of Scenes: ", len(scenes))

# The first scene has all the bands
scene = scenes[0]

scene.load(['C10'])
new_scn = scene.resample(scene.min_area(), resampler='native')

with ProgressBar():
    new_scn.save_datasets(filename='test.tif', overlay={'coast_dir': '/tmp/shapefiles/', 'color': (0, 0, 0), 'width': 3, 'level_coast': 1, 'level_borders': 2})

Note that the same code without the overlay argument in save_datasets works without issue

And can you do print(scn['C10'])?

<xarray.DataArray 'Rad' (y: 5424, x: 5424)>
dask.array<shape=(5424, 5424), dtype=float64, chunksize=(4096, 4096)>
Coordinates:
  * x        (x) float64 -0.1518 -0.1518 -0.1517 ... 0.1517 0.1518 0.1518
    y_image  float32 0.0
    x_image  float32 0.0
  * y        (y) float64 0.1518 0.1518 0.1517 0.1517 ... -0.1517 -0.1518 -0.1518
Attributes:
    satellite_longitude:    -75.19999694824219
    satellite_latitude:     0.0
    satellite_altitude:     35786.0234375
    standard_name:          toa_brightness_temperature
    _FillValue:             4095
    _Unsigned:              true
    modifiers:              ()
    grid_mapping:           goes_imager_projection
    add_offset:             -0.9561
    name:                   C10
    sensor:                 abi
    ancillary_variables:    []
    wavelength:             (7.24, 7.34, 7.44)
    platform_name:          GOES-16
    scale_factor:           0.02004128
    valid_range:            [   0 4094]
    units:                  K
    long_name:              ABI L1b Radiances
    resolution:             2000
    calibration:            brightness_temperature
    sensor_band_bit_depth:  12
    cell_methods:           t: point area: point
    start_time:             2018-10-07 16:00:38
    end_time:               2018-10-07 16:11:15.900000
    area:                   Area ID: abi_geos\nDescription: ABI L1B file area...
    polarization:           None
    level:                  None

Side question: How do you download those files from AWS without an AWS client?

I used the aws-cli to download each of the nc files.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Hm I wonder if the native resampling is messing with this. FYI for a single channel you don't need to use the native resampler. That is only for composites trying to use data in different resolutions.

Edit: Nope, nevermind. I was able to reproduce.

@hany
Copy link
Author

hany commented Oct 16, 2018

Good to know. Just tried without resampling, same error:

Traceback (most recent call last):
  File "testsatpy.py", line 36, in <module>
    scene.save_datasets(overlay={'coast_dir': '/tmp/shapefiles/', 'color': (0, 0, 0), 'width': 3, 'level_coast': 1, 'level_borders': 2})
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/scene.py", line 1071, in save_datasets
    return writer.save_datasets(datasets, compute=compute, **save_kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 462, in save_datasets
    res = self.save_dataset(ds, compute=False, **kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 565, in save_dataset
    fill_value=fill_value, **kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/geotiff.py", line 192, in save_image
    filename = filename or self.get_filename(**img.data.attrs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/satpy/writers/__init__.py", line 427, in get_filename
    return self.filename_parser.compose(kwargs)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/trollsift/parser.py", line 55, in compose
    return compose(self.fmt, keyvals)
  File "/Users/hany/.virtualenvs/goes16-py3/lib/python3.6/site-packages/trollsift/parser.py", line 355, in compose
    return formatter.format(fmt, **keyvals)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 190, in format
    return self.vformat(format_string, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 194, in vformat
    result, _ = self._vformat(format_string, args, kwargs, used_args, 2)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 234, in _vformat
    obj, arg_used = self.get_field(field_name, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 299, in get_field
    obj = self.get_value(first, args, kwargs)
  File "/Users/hany/.pyenv/versions/3.6.5/lib/python3.6/string.py", line 256, in get_value
    return kwargs[key]
KeyError: 'name'

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Note that right now the default enhancements and limitations of pycoast will not let you add coastlines to a single band dataset. See #449 for details and status. So even after I fix this KeyError the single band will still fail.

@hany
Copy link
Author

hany commented Oct 16, 2018

Strange, it actually works for me (assuming I specify the filename argument to save_datasets). I do not get the error from #449.

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Interesting...very interesting...I'll play around with it and see what I get. What version of pycoast, trollimage, and satpy are you using?

@hany
Copy link
Author

hany commented Oct 16, 2018

All installed via pip:

pycoast==1.1.0
satpy==0.9.4
trollimage==1.5.7

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

What version of aggdraw?

@hany
Copy link
Author

hany commented Oct 16, 2018

Also installed via pip:

aggdraw==1.3.7

I've attached my requirements.txt in case it helps.
requirements.txt

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Oh you are colorizing your datasets aren't you? You have a custom enhancement, right?

@hany
Copy link
Author

hany commented Oct 16, 2018

Correct, from the efforts in #459. Let me try without.

@hany
Copy link
Author

hany commented Oct 16, 2018

OK, without colorizing, I do get the error from #449!

@djhoese
Copy link
Member

djhoese commented Oct 16, 2018

Ok this all makes sense now. The issue is that pycoast doesn't know what to do with single band (mode='L') data because it assumes that colors provided by the user are RGB. So as mentioned in #449 you either have to convert to RGB before hand or fix pycoast to support this.

Since you are colorizing your data they are being converted from L to RGB so you were getting around the issue.

@hany
Copy link
Author

hany commented Oct 16, 2018

Makes perfect sense!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants