Pycoast v0.5.2 issue when drawing coastlines

[sassimarco]

Reference to similar problem in Google Group:
https://groups.google.com/forum/#!searchin/pytroll/pycoast$20lines%7Csort:date/pytroll/tzmtaDuaGD4/kVrY2RZ93gEJ

Hrobjartur Thorsteinsson:
Coastlines need a consistency check that drawn line-segments do not exceed some sort of threshold. We had considered that, but now I think I should add that improvement.

[djhoese]

I know this is old, but what is the actual version pycoast version at fault here? The current version released is 1.2.1.

[djhoese]

Ah it is 0.5.2 as mentioned in the linked google thread.

[djhoese]

Could someone (@sassimarco, @peters77, ?) make a self-contained example for this one? I just merged some stuff to master that handles thresholds better, but this was really only intended for some pyproj changes. Probably not going to fix this, but hard to tell without an example that doesn't require outside data.

[sassimarco]

Hi Dave, here an example, quite an old implementation.... [image: image.png] dask==1.0.0 h5netcdf==0.5.0 h5py==2.8.0 netCDF4==1.4.2 numpy==1.15.2 pandas==0.23.4 path.py==11.0 Pillow==5.2.0 pycoast==1.1.0 pykdtree==1.3.1 pyninjotiff==0.1.0 pyproj==1.9.5.1 pyresample==1.10.3 pyshp==1.2.12 pyspectral==0.8.3 python-dateutil==2.7.3 python-geotiepoints==1.1.6 pytz==2018.5 PyYAML==3.13 requests==2.19.1 satpy==0+unknown scipy==1.1.0 six==1.11.0 toolz==0.9.0 trollimage==1.5.7 trollsift==0.2.1 urllib3==1.23 xarray==0.11.0

…

On Fri, 6 Dec 2019 at 21:42, David Hoese ***@***.***> wrote: Could someone ***@***.*** </~https://github.com/sassimarco>, @peters77 </~https://github.com/peters77>, ?) make a self-contained example for this one? I just merged some stuff to master that handles thresholds better, but this was really only intended for some pyproj changes. Probably not going to fix this, but hard to tell without an example that doesn't require outside data. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#23?email_source=notifications&email_token=AK2I2LZNNARGVLZR5IS5QYTQXK2KHA5CNFSM4GGUMDXKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEGFJJBY#issuecomment-562730119>, or unsubscribe </~https://github.com/notifications/unsubscribe-auth/AK2I2LYSQ2U6DHTQ62SAUC3QXK2KHANCNFSM4GGUMDXA> .

[djhoese]

@sassimarco I don't think github liked you sending an image over email. You may have to add the image directly from the github page. If you have the code as text that I can copy and paste that would be the best kind of example.

[sassimarco]

[sassimarco]

#!/usr/bin/python
#
#  python  geos17_processing.py

# general modules
# -------------------------
import glob
import os
import sys
import time, datetime
from path import Path
from PIL import Image, ImageFont, ImageDraw, PngImagePlugin
import numpy as np
import shutil

# Pytroll modules
# -------------------------
from satpy.scene import Scene
from satpy.utils import debug_on
from satpy.writers import to_image, get_enhanced_image
import pyninjotiff
from pycoast import ContourWriter

# debug mode for Pytroll modules
# -------------------------
debug_on()
import warnings
warnings.filterwarnings("ignore")

# defines working variables
# -------------------------
goes_input_dir      = '/var/tmp/goes17/'
goes_PNG_output_dir = '/opt/pytroll/outbox/intranet/'
goes_TIF_output_dir = '/opt/pytroll/outbox/ninjo/'
goes_TIF_PROD_dir   = '/opt/pytroll/outbox/ninjo_prod/'
shapes_dir          = '/opt/pytroll/shapes/'
tiffdump_file       = '/var/tmp/goes17.tif'

sat_live = True

print("")
print("---------------------------------------")
print("---------- GOES-17 processing ---------")
print("---------------------------------------")
print("")

# extract timestamp from last incoming ABI file
# ---------------------------------------------
ABI_files = glob.glob(goes_input_dir + 'OR_ABI-L1b*.nc')

if not ABI_files:
    print("===============================================")
    print("GOES-17 data not available, skip processing....")
    print("===============================================")
    sys.exit()

latest_ABI_file = os.path.basename(max(ABI_files, key=os.path.getctime))
print("latest GOES17 file-> ", latest_ABI_file)

starttime = (latest_ABI_file[latest_ABI_file.find("s")+1:latest_ABI_file.find("_e")])
year   = starttime[0:4]
julian = starttime[4:7]
hours  = starttime[7:9]
min    = starttime[9:11]
stime  = datetime.datetime(int(year), 1, 1) + datetime.timedelta(int(julian)-1) + datetime.timedelta(hours=int(hours)) + datetime.timedelta(minutes=int(min))
print("starttime --> ", stime)

# timestamp for final products
timestamp = year + str(stime.month).zfill(2) + str(stime.day).zfill(2) + hours + min

endtime = (latest_ABI_file[latest_ABI_file.find("e")+1:latest_ABI_file.find("_c")])
year   = endtime[0:4]
julian = endtime[4:7]
hours  = endtime[7:9]
min    = endtime[9:11]
etime  = datetime.datetime(int(year), 1, 1) + datetime.timedelta(int(julian)-1) + datetime.timedelta(hours=int(hours)) + datetime.timedelta(minutes=int(min))
print("endtime     --> ", etime)

# select channels
# -----------------------------

channels = ['C02', 'C08', 'C13']
channel_list_to_process = [] 

channels_to_ninjo = {}
channels_to_ninjo['C02'] = [ 'GOES-17 channel C02', 'VIS006', 'nrGOESW8km']
channels_to_ninjo['C08'] = [ 'GOES-17 channel C08', 'WV062', 'nrGOESW8km']
channels_to_ninjo['C13'] = [ 'GOES-17 channel C13', 'IR107', 'nrGOESW8km']

# check which channel is available, with a correct timestamp
for goes17_channel in channels:
    pattern = 'OR_ABI-L1b*M?' + goes17_channel + '*' + starttime + '*.nc'
    if glob.glob(goes_input_dir + pattern):
        print ("GOES-17 data channel " + goes17_channel + " is available")
        channel_list_to_process.append(goes17_channel)
    else:
        print ("GOES-17 data channel " + goes17_channel + " is not available")


# exit if there are no products to process
# ---------------------------------------
if not channel_list_to_process:
    print("=======================================")
    print("no realtime GOES-17 data to process")
    sys.exit("no products to process --> exit")

# exit if data timestamp older is than 1 hour
# -------------------------------------------
now = datetime.datetime.now()
if now - datetime.timedelta(minutes=60) <= stime <= now:
    print "GOES-17 data are on time"
else:
    print("===================================================")
    sys.exit("GOES-17 data are older than 1 hour --> exit")


# build scene 
# -------------------------------
try:
#    scn = Scene(
#        platform_name = "GOES-17",   
#        sensor = "abi",
#        start_time = stime,
#        end_time = etime,
#        base_dir = goes_input_dir,
#        reader = 'abi_l1b'
#    )

    scn = Scene(reader="abi_l1b", filenames=ABI_files)

except:
    print("===================================================")
    sys.exit("GOES-17 data could not be loaded (different timestamp?)")

cw = ContourWriter(shapes_dir)

# start resampling and produce output  
# ------------------------------------
from satpy.resample import get_area_def
areaid = 'NinJoGOESWregion'
areadef = get_area_def(areaid)

for channel in channel_list_to_process:
    print ("Processing channel " + channel + "...")
    # file name build on Cinesat  "pattern" --> GOES15-W_IR107_nrGOESW8km_1709130700.tif
    outfile = 'GOES17-W_' + channels_to_ninjo.get(channel)[1] + '_' + channels_to_ninjo.get(channel)[2] + '_' + timestamp[2:]
    
    try:
        scn.load([channel])
    except:
        print("===================================================")
        sys.exit("GOES-17 data " + channel + " could not be loaded (different timestamp?)")


    lcd = scn.resample(areadef, cache_dir='/var/tmp/sam/')

    if channel == 'C02':
        lcd.save_dataset(channel, filename = tiffdump_file,
                         writer='ninjotiff',
                         ninjo_product_name = 'G17_'+ channel,
                         zero_seconds=True)
    else:
        # Min and max values in Celsius
        value_min = -87.5
        value_max = 40

     #   lcd[channel].clip(value_min - 273.15, value_max - 273.15)

        # physic_unit='C' is used to invert the color table
        lcd.save_dataset(channel, filename = tiffdump_file,
                         writer='ninjotiff',
                         ninjo_product_name = 'G17_'+ channel,
                         physic_unit='K',
                         zero_seconds=True,
                         enhancement_config=False,
                         fill_value=255,
                         ch_max_measurement_unit=value_max,
                         ch_min_measurement_unit=value_min)
                         
    
    # make a copy just for NinJO....
    shutil.copyfile(tiffdump_file, goes_TIF_output_dir + outfile + '.tif')
    print('copied to --> ', goes_TIF_output_dir + outfile + '.tif')

    # make a copy for the PROD delivery to NinJo
    shutil.copyfile(tiffdump_file, goes_TIF_PROD_dir + outfile + '.tif')
    print('copied to --> ', goes_TIF_PROD_dir + outfile + '.tif')

    # save as PNG graphic file
    img = get_enhanced_image(lcd[channel])
    if channel <> 'C02': img.invert()
    # stretching needed because was deactivated as "dummy" in satpy/etc/enhancements/generic.yaml to avoid "stretching" for ninjotiff
    img.save(goes_PNG_output_dir + outfile + '.png')

    # add borders & coasts
    img = Image.open(goes_PNG_output_dir + outfile  + '.png')
    cw.add_coastlines(img, areadef)
    cw.add_borders(img, areadef, resolution='i')

    # save as png graphic file
    img.save(goes_PNG_output_dir + outfile + '.png')

    # GOES-17 images resized only for the satellite LIVE
    # --------------------------------------------------
    if sat_live:
        basewidth=1200
        wpercent = (basewidth/float(img.size[0]))
        hsize = int((float(img.size[1])*float(wpercent)))
        img = img.resize((basewidth, hsize), Image.ANTIALIAS)
        img = img.crop((0+basewidth/7, 0+hsize/7, basewidth-basewidth/7, hsize-hsize/7))
        img.convert('RGB').save(goes_PNG_output_dir + outfile + '.jpg')

# ----------------------
# cleanup
# -----------------------
d = Path(goes_PNG_output_dir)

removed = 0
hours = 6
time_in_secs = time.time() - hours*3600

for file in d.walk():
    if file.isfile():
        if file.mtime <= time_in_secs:
            file.remove()
            removed += 1

print "---------"
print "removed old GOES-17 products --> ", removed
print("---------------------------------------")

[djhoese]

I was kind of hoping for a simple example (15 lines?). I'm hoping that the data can be randomly generated with np.random.random and the AreaDefinition will be hardcoded.

[sassimarco]

in that case these are the lines...

    # add borders & coasts
    img = Image.open(goes_PNG_output_dir + outfile  + '.png')
    cw.add_coastlines(img, areadef)
    cw.add_borders(img, areadef, resolution='i')

    # save as png graphic file
    img.save(goes_PNG_output_dir + outfile + '.png')

[djhoese]

Thanks!

FYI I updated your last two comments. The code blocks need to have three backticks before and after the code. The three backticks have to be on their own lines though (no code on those lines).

[mraspaud]

@sassimarco what's the area definition you use ? Can you paste it here ?

[sassimarco]

Ciao Martin, ok, i found it out in a old GOES17 Script..... areaid = 'NinJoGOESWregion' areadef = get_area_def(areaid) REGION: NinJoGOESEregion { NAME: NinJoGOESEregion PCS_ID: NinJoGOESEregion PCS_DEF: proj=eqc, lon_0=-75.2, ellps=WGS84 XSIZE: 2500 YSIZE: 2500 AREA_EXTENT: (-10008769.814946,-10016778.695711,10016778.695711,10008768.924390) }; REGION: NinJoGOESWregion { NAME: NinJoGOESWregion PCS_ID: NinJoGOESWregion PCS_DEF: proj=eqc, lon_0=-137.2, a=6378137.00 XSIZE: 2500 YSIZE: 2500 AREA_EXTENT: (-10008769.814946,-10016778.695711,10016778.695711,10008768.924390) }; Marco

…

On Fri, 21 Feb 2020 at 12:12, Martin Raspaud ***@***.***> wrote: @sassimarco </~https://github.com/sassimarco> what's the area definition you use ? Can you paste it here ? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#23?email_source=notifications&email_token=AK2I2L7NN7WVIVYIGWNQJITRD6ZLHA5CNFSM4GGUMDXKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMSLYQY#issuecomment-589610051>, or unsubscribe </~https://github.com/notifications/unsubscribe-auth/AK2I2L25EUGX5MVGZKVCVQTRD6ZLHANCNFSM4GGUMDXA> .

[mraspaud]

Thanks.

[gerritholl]

Related: pytroll/satpy#1224 contains a MCVE that with pytroll/satpy@aed86fd reproduces the problem (a roundabout way, but still) when reading an ICON GRIB file:

from satpy import Scene
from glob import glob
sc = Scene(filenames={"grib":
    glob("/media/nas/x21308/scratch/ICON/*12:00:00Z*000.grib")})
sc.load(["t"])
sc.show("t", overlay={"coast_dir": "/media/nas/x21308/shp/", "color": "red"})

which with the linked pull request gives

[yzhan-met]

Hi there, I tried to draw coastline for the given area definition and found this issue as well. Is there an existing way to get around it? Thanks

Here are my area definition and script.

areas.yaml

TROPICS:
  projection:
    proj: merc
    lat_0: -13.5
    lon_0: 177.5
    ellps: WGS84
  resolution: [2000, 2000]
  area_extent:
    upper_right_xy: [4175000, 1097830]
    lower_left_xy: [-4175000, -4140170]
    units: m

In [43]: import pyresample
In [44]: from PIL import Image
In [45]: from pycoast import ContourWriterAGG
In [46]: cw = ContourWriterAGG('/tmp/gshhs')
In [47]: img = Image.new('RGB', (425, 425))
    ...: area_def = pyresample.load_area("/tmp/areas.yaml", "Tropics")
    ...: cw.add_coastlines(img, area_def, resolution='i')
    ...: img.show()

[djhoese]

Thanks for the example @yzhan-met. I simplified it a bit with modern pyresample so it doesn't require the extra YAML:

import pyresample
from PIL import Image
from pycoast import ContourWriterAGG
cw = ContourWriterAGG('/data/gshhg_shapefiles')
img = Image.new('RGB', (425, 425))
area_def = pyresample.create_area_def(
    "test",
    "+proj=merc +lat_0=-13.5 +lon_0=177.5 +ellps=WGS84",
    width=425, height=425,
    area_extent=[-4175000, -4140170, 4175000, 1097830],
)
cw.add_coastlines(img, area_def, resolution='i')
img.show()

We've had a user on the mailing list come up with a workaround:

--- cw_base.py.dist	2021-08-18 17:14:34.000000000 +0000
+++ cw_base.py	2022-03-11 09:41:04.664157608 +0000
@@ -593,9 +593,14 @@
 
             # Check if polygon is possibly relevant
             s_lon_ll, s_lat_ll, s_lon_ur, s_lat_ur = shape.bbox
-            shape_is_outside_lon = lon_max < s_lon_ll or lon_min > s_lon_ur
+            if lon_min <= lon_max:
+                # Area_extent west or east of dateline
+                shape_is_outside_lon = lon_max < s_lon_ll or lon_min > s_lon_ur
+            else:
+                # Area_extent spans over dateline
+                shape_is_outside_lon = lon_max < s_lon_ll and lon_min > s_lon_ur
             shape_is_outside_lat = lat_max < s_lat_ll or lat_min > s_lat_ur
-            if lon_min <= lon_max and (shape_is_outside_lon or shape_is_outside_lat):
+            if shape_is_outside_lon or shape_is_outside_lat:
                 # Polygon is irrelevant
                 continue

This doesn't solve every case of lines going across the image but it certainly helps. I just wanted to record this somewhere other than the mailing list in case people had other ideas or other solutions. I also plan on bringing up other ideas with the user on the mailing list.