readshapefile() inconsistencies

[avipersin]

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='robin', lat_0=0, lon_0=100)
map.readshapefile('./Paleo_Sturtian_750Ma', 'Paleo_Sturtian_750Mm')
plt.show()

The shapefile I'm using can be downloaded from: my github.

Using the "robin" projection:

• lon_0=0 everything works fine
• lon_0 = anything else; causes wrap around problems. Horizontal lines through the plot.
  

Using the "cyl" projection:

• lon_0=0 everything works fine.
• lon_0 = anything else; No lines through the plot but doesn't plot the right side of the shape file.
  

Is it possible to change the center longitude while plotting a shapefile?

[WeatherGod]

The problem is that it is difficult to unwrap the coordinates for polygons because it depends upon which way across the fold the polygon is going. This is one of the reasons why basemap is being replaced by cartopy, which handles situations like these better.

…

On Wed, Jun 13, 2018 at 1:33 PM, avipersin ***@***.***> wrote: from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt map = Basemap(projection='robin', lat_0=0, lon_0=100) map.readshapefile('./Paleo_Sturtian_750Ma', 'Paleo_Sturtian_750Mm') plt.show() Using the "robin" projection: - lon_0=0 everything works fine - lon_0 = anything else; causes wrap around problems. Horizontal lines through the plot. [image: robin] <https://user-images.githubusercontent.com/8727319/41367816-224d3206-6f0e-11e8-9a17-7b994f972bf5.png> Using the "cyl" projection: - lon_0=0 everything works fine. - lon_0 = anything else; No lines through the plot but doesn't plot the right side of the shape file. [image: cyl] <https://user-images.githubusercontent.com/8727319/41367795-12d6974a-6f0e-11e8-8d1b-35612f363983.png> Is it possible to change the center longitude while plotting a shapefile? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#410>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARy-CsIyx29NRy0OlEjd-nVzmeMDMvCks5t8UzTgaJpZM4UmoLC> .

[avipersin]

Is there a workaround or a continental outline format supported by basemap that doesn't have this issue?

How does drawcoastlines() avoid this problem entirely?

[WeatherGod]

to be honest, I am actually not certain that it does properly avoid the situation. We do have a shiftdata function, but it is finicky at best and can produce wrong results. I don't know for sure if drawcoastlines() uses it or not.

[avipersin]

In case others stumble upon this...

The readshapefile() function is glitchy. There are inconsistencies when using different projections as well as different center longitude and latitude points. In order to get this working I've had to convert my original shapefiles to one that the readshapefile() function can properly understand.

The code to convert the shapefiles:

import matplotlib as mpl
from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import addcyclic
import shapefile
import math

path_to_shapefile = "path_to_shapefile"
cen_lon = -100
cen_lat = -45
projection = "ortho"
bm = Basemap(projection=projection, lon_0=cen_lon, lat_0=cen_lat)

w = shapefile.Writer()
w.field('name', 'C')
sf = shapefile.Reader(path_to_shape_file)
for shpid, shape in enumerate(sf.shapes()):
    if shpid > -1:
        points = shape.points
        line = []
        lines = []
        for i in range(len(points) - 1):
            # if i < 21:
            #     continue
            lon1 = points[i][0]
            lon2 = points[i + 1][0]
            lat1 = points[i][1]
            lat2 = points[i + 1][1]

            if lon1 < (cen_lon - 180):
                lon1 = lon1 + 360
            elif lon1 > (cen_lon + 180):
                lon1 = lon1 - 360

            if lon2 < (cen_lon - 180):
                lon2 = lon2 + 360
            elif lon2 > (cen_lon + 180):
                lon2 = lon2 - 360

            if lon2 < (cen_lon + 180) < lon1:
                if lon1 + lon2 < 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon - 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon + 180, (lat1 + lat2) / 2)]
                elif lon1 + lon2 >= 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon + 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon - 180, (lat1 + lat2) / 2)]
            elif lon1 < (cen_lon + 180) < lon2:
                if lon1 + lon2 < 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon - 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon + 180, (lat1 + lat2) / 2)]
                elif lon1 + lon2 >= 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon + 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon - 180, (lat1 + lat2) / 2)]
            elif lon1 <= (cen_lon + 180) and abs(lon1 - lon2) > 350:
                if lon1 - lon2 >= 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon + 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon - 180, (lat1 + lat2) / 2)]
                elif lon1 - lon2 < 0:
                    line.append((lon1, lat1))
                    line.append((cen_lon - 180, (lat1 + lat2) / 2))
                    lines.append(line)
                    line = [(cen_lon + 180, (lat1 + lat2) / 2)]
            else:
                line.append((lon1, lat1))

        if shape.points[-1][0] < (cen_lon - 180):
            lon1 = shape.points[-1][0] + 360
        elif shape.points[-1][0] > (cen_lon + 180):
            lon1 = shape.points[-1][0] - 360
        else:
            lon1 = shape.points[-1][0]
        line.append((lon1, shape.points[-1][1]))
        lines.append(line)

        if projection in ["ortho"]:
            line_new = []
            for line in lines:
                for i in range(len(line)):
                    lonC = cen_lon
                    latC = cen_lat
                    phiCRad = latC * math.pi / 180.
                    cosPhiC = math.cos(phiCRad)
                    sinPhiC = math.sin(phiCRad)
                    lon = line[i][0]
                    lat = line[i][1]
                    dlon = lon - lonC
                    if dlon > 180:
                        dlon -= 360
                    elif dlon < -180:
                        dlon += 360
                    lambdaRad = dlon * math.pi / 180
                    cosLambda = math.cos(lambdaRad)

                    phiRad = lat * math.pi / 180
                    cosPhi = math.cos(phiRad)
                    sinPhi = math.sin(phiRad)

                    cosZ = sinPhiC * sinPhi + cosPhiC * cosPhi * cosLambda

                    if cosZ < 0:
                        continue
                    elif round(cosZ, 5) == 0:
                        continue
                    else:
                        line_new.append((lon, lat))

            tmp_lines = []
            tmp_i = 0
            if len(line_new) > 0:
                for i in range(len(line_new) - 1):
                    if abs(line_new[i][0] - line_new[i + 1][0]) >= 10 or abs(
                            line_new[i][1] - line_new[i + 1][1]) >= 10:
                        tmp_lines.append(line_new[tmp_i:i + 1])
                        tmp_i = i + 1
                tmp_lines.append(line_new[tmp_i:])

                for lines in tmp_lines:
                    w.record('line' + str(shpid))
                    w.line([lines])

        else:
            for idx, shape in enumerate(lines):
                w.record('line' + str(shpid) + "_" + str(idx))
                w.line([shape])
        w.save('tmp_shapefile')

bm.readshapefile('tmp_shapefile', 'tmp_shapefile')

plt.show()