Soru matplotlib / Basemap ile nehirler olmadan dünya haritası?


Kıtaların sınırlarını Basemap'la (ya da başka bir yol varsa, Basemap olmadan), bu can sıkıcı nehirler gelmeden çizmenin bir yolu var mıdır? Özellikle de Kongo Nehrinin, hatta okyanusa bile ulaşmadığı bir parça rahatsız ediyor.

DÜZENLEME: Haritadaki verileri Basemap galerisi (ve yine de, verilerin üzerinde siyah çizgiler olarak çizilmiş kıtaların sınır çizgileri var, dünya haritasına bir yapı kazandırmak için) bu yüzden Kanca ile çözüm, güzel, ustaca bile olsa, bu amaç için geçerli değildir.

world map

Görüntü tarafından üretilen:

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

fig = plt.figure(figsize=(8, 4.5))
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.00)
m = Basemap(projection='robin',lon_0=0,resolution='c')
m.fillcontinents(color='gray',lake_color='white')
m.drawcoastlines()
plt.savefig('world.png',dpi=75)

18
2018-01-11 14:34


Menşei


İltifatınız için çok teşekkürler! İhtiyaçlarınıza en uygun cevabı seçmek için çekinmeyin (bu site ne işe yarar!). Gelecekte başka birine yararlı olabileceğinden, cevabımı bırakacağım. Gelecekte belirtmeye çalışın kesinlikle ne yapmak istersin Sunulan iki cevap, her ikisini de sanatsal bir yorum olarak yorumladığımız "komplo" bölümünü ele alıyor. Bununla birlikte, ikimizin de sunduğumuz cevabı hala çizemediğin için bir neden yok - tüm koordinatlar orada, şu anda bütün bu çirkin nehirlerin nerede olduğunu biliyorsun! - Hooked


Cevaplar:


Bunun gibi nedenlerden ötürü sıklıkla Basemap'tan kaçınıyorum ve OGR ile şekil dosyasını okuyup bir Matplotlib sanatçısına dönüştürüyorum. Hangisi daha çok iş ama aynı zamanda daha fazla esneklik sağlıyor.

Temel harita, girdi verilerinin koordinatlarını 'çalışma projeksiyonunuza' dönüştürmek gibi bazı çok iyi özelliklere sahiptir.

Basemap ile sopa yapmak istiyorsanız, nehirleri içermeyen bir şekil dosyası alın. Örneğin Doğal Dünya, fiziksel bölümde ('ölçek rank' verisini indir ve aç) 'güzel bir' Land 'şekil dosyasına sahiptir. Görmek http://www.naturalearthdata.com/downloads/10m-physical-vectors/

Şekil dosyasını Basemap'ten m.readshapefile () yöntemi ile okuyabilirsiniz. Bu, Matplotlib Path köşelerini ve kodlarını, daha sonra yeni bir Yol'a dönüştürebileceğiniz projeksiyon koordinatlarına almanızı sağlar. Onun biraz dolambaçlı ama Matplotlib'den, çoğu Basemap üzerinden doğrudan erişilemeyen tüm stil seçenekleri sunuyor. Onun biraz hackish, ama Basemap'a yapışırken başka bir şekilde değilim.

Yani:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.collections import PathCollection
from matplotlib.path import Path

fig = plt.figure(figsize=(8, 4.5))
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.00)

# MPL searches for ne_10m_land.shp in the directory 'D:\\ne_10m_land'
m = Basemap(projection='robin',lon_0=0,resolution='c')
shp_info = m.readshapefile('D:\\ne_10m_land', 'scalerank', drawbounds=True)
ax = plt.gca()
ax.cla()

paths = []
for line in shp_info[4]._paths:
    paths.append(Path(line.vertices, codes=line.codes))

coll = PathCollection(paths, linewidths=0, facecolors='grey', zorder=2)

m = Basemap(projection='robin',lon_0=0,resolution='c')
# drawing something seems necessary to 'initiate' the map properly
m.drawcoastlines(color='white', zorder=0)

ax = plt.gca()
ax.add_collection(coll)

plt.savefig('world.png',dpi=75)

verir:

enter image description here


9
2018-01-11 16:11



+1 Doğal Dünya referansı için teşekkür ederim! "Sanatsal" render için başka iyi açık kaynaklar biliyor musunuz? - Hooked
Diğer Matplotlib bazen ben Mapnik (TileMill ile) kullanıyorum. Ben sadece yukarıdaki yöntem benim çokgenler iç halkalar ihmal / kaldırır fark ettim (yani Kaspian denizi örneğin dissapears) :( - Rutger Kassies
Şimdiye kadar teşekkürler, ama: ne_10m_land.zip linkinden indirdim, unzipped ve 'ne_10m_land' olarak okudum, ama hatalar alıyorum: "satır 10, <module> shp_info = m.readshapefile ('ne_10m_land', 'scalerank', drawbounds = True) [...] ValueError: base 10 ile int () için geçersiz literal: '****' ". - Sampo Smolander
Bence onun orijinal şekil dosyası 'scalerank' özniteliğinde geçerli tamsayılar içermiyor. Bazı sahte numaralar girmeye yardımcı olabilir. Herhangi bir özellik yoksa Basemap'in 'fid' seçip seçemeyeceğini bilmiyorum, güzel olurdu. - Rutger Kassies
Bir temel harita alternatifi tavsiye edebilir misiniz? Abandonware gibi görünüyor. Teşekkürler. - tommy.carstensen


"Can sıkıcı" nehirler nasıl kaldırılır:

Resmi post-proseslemek istiyorsanız (doğrudan Basemap ile çalışmak yerine) okyanusa bağlanmayan su kütlelerini kaldırabilirsiniz:

import pylab as plt
A = plt.imread("world.png")

import numpy as np
import scipy.ndimage as nd
import collections

# Get a counter of the greyscale colors
a      = A[:,:,0]
colors = collections.Counter(a.ravel())
outside_and_water_color, land_color = colors.most_common(2)

# Find the contigous landmass
land_idx = a == land_color[0]

# Index these land masses
L = np.zeros(a.shape,dtype=int) 
L[land_idx] = 1
L,mass_count = nd.measurements.label(L)

# Loop over the land masses and fill the "holes"
# (rivers without outlays)
L2 = np.zeros(a.shape,dtype=int) 
L2[land_idx] = 1
L2 = nd.morphology.binary_fill_holes(L2)

# Remap onto original image
new_land = L2==1
A2 = A.copy()
c = [land_color[0],]*3 + [1,]
A2[new_land] = land_color[0]

# Plot results
plt.subplot(221)
plt.imshow(A)
plt.axis('off')

plt.subplot(222)
plt.axis('off')
B = A.copy()
B[land_idx] = [1,0,0,1]
plt.imshow(B)

plt.subplot(223)
L = L.astype(float)
L[L==0] = None
plt.axis('off')
plt.imshow(L)

plt.subplot(224)
plt.axis('off')
plt.imshow(A2)

plt.tight_layout()  # Only with newer matplotlib
plt.show()

enter image description here

İlk görüntü orjinaldir, ikincisi kara kütlesini tanımlar. Üçüncüsü gerekli değil ama ID'nin her bir ayrı bitişik toprak kütlesi kadar eğlenceli. Dördüncü resim istediğin şey, "nehirler" olan görüntü kaldırıldı.


6
2018-01-11 16:10



Dördüncü fotoğrafta hala nehirler var. - tiago
@tiago Doğru, ama cevabımı okuduysanız, sunulan yöntemin "okyanusun akıntısı olmayan su kütlelerini" çıkarmak olduğunu fark etmiş olursunuz, ki bu OP'nin "can sıkıcı nehirler" olarak adlandırdığı farzedildi. Gözlemler yapmak açıkça belirtilmiş söylenen kadar yapıcı değildir, kalan eksiklikleri düzeltmek için bir öneri. - Hooked


User1868739 örneğini takip ederek, yalnızca istediğim yolları (bazı göller için) seçebiliyorum: world2

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

fig = plt.figure(figsize=(8, 4.5))
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.00)
m = Basemap(resolution='c',projection='robin',lon_0=0)
m.fillcontinents(color='white',lake_color='white',zorder=2)
coasts = m.drawcoastlines(zorder=1,color='white',linewidth=0)
coasts_paths = coasts.get_paths()

ipolygons = range(83) + [84] # want Baikal, but not Tanganyika
# 80 = Superior+Michigan+Huron, 81 = Victoria, 82 = Aral, 83 = Tanganyika,
# 84 = Baikal, 85 = Great Bear, 86 = Great Slave, 87 = Nyasa, 88 = Erie
# 89 = Winnipeg, 90 = Ontario
for ipoly in ipolygons:
    r = coasts_paths[ipoly]
    # Convert into lon/lat vertices
    polygon_vertices = [(vertex[0],vertex[1]) for (vertex,code) in
                        r.iter_segments(simplify=False)]
    px = [polygon_vertices[i][0] for i in xrange(len(polygon_vertices))]
    py = [polygon_vertices[i][2] for i in xrange(len(polygon_vertices))]
    m.plot(px,py,linewidth=0.5,zorder=3,color='black')

plt.savefig('world2.png',dpi=100)

Ancak bu sadece kıtalar için beyaz arka plan kullanırken çalışır. Değiştirirsem color için 'gray' Aşağıdaki satırda, diğer nehirlerin ve göllerin kıtalarla aynı renkle doldurulmadığını görüyoruz. (Ayrıca area_thresh okyanusa bağlı olan nehirleri çıkarmayacak.

m.fillcontinents(color='gray',lake_color='white',zorder=2)

world3

Beyaz arka plana sahip versiyon, kıtalar üzerindeki her türlü arazi bilgisini renklendirmek için yeterlidir, ancak kıtalar için gri arkaplanı korumak isteyen daha ayrıntılı bir çözüme ihtiyaç olacaktır.


3
2018-01-18 14:11



Eğer "istenmeyen" çokgenler çizerseniz ve aynı zamanda kıta dolgusuyla aynı renkte olursanız beyaz ırmak probleminin ortadan kalkacağını tahmin ediyorum. - user1868739


Bu 'kırık' nehirlerden kaçınmak için Basemap'ın drawcoastlines () 'ını sık sık değiştiriyorum. Ayrıca veri kaynağı tutarlılığı için drawcountries () değiştiriyorum.

Natural Earth verilerinde bulunan farklı çözünürlükleri desteklemek için kullandığım şey şöyledir:

from mpl_toolkits.basemap import Basemap


class Basemap(Basemap):
    """ Modify Basemap to use Natural Earth data instead of GSHHG data """
    def drawcoastlines(self):
        shapefile = 'data/naturalearth/coastline/ne_%sm_coastline' % \
                    {'l':110, 'm':50, 'h':10}[self.resolution]
        self.readshapefile(shapefile, 'coastline', linewidth=1.)
    def drawcountries(self):
        shapefile = 'data/naturalearth/countries/ne_%sm_admin_0_countries' % \
                    {'l':110, 'm':50, 'h':10}[self.resolution]
        self.readshapefile(shapefile, 'countries', linewidth=0.5)


m = Basemap(llcrnrlon=-90, llcrnrlat=-40, urcrnrlon=-30, urcrnrlat=+20,
            resolution='l')  # resolution = (l)ow | (m)edium | (h)igh
m.drawcoastlines()
m.drawcountries()

İşte çıktı

Lütfen, varsayılan olarak Basemap'in gösterilen kodda desteklenmeyen çözünürlük = 'c' (ham) olduğunu unutmayın.


3
2018-03-10 19:20



Bu çözüm için verilen verileri indirmeniz gerekir. Doğal toprak. - safay


Şekil dosyalarından ziyade anahatları çizmek için sorun yoksa, her yerden alabileceğiniz sahil hatlarını çizmek oldukça kolaydır. Benim sahil hatlarını MATLAB formatında NOAA Sahil şeridi Extractor'tan aldım:      http://www.ngdc.noaa.gov/mgg/shorelines/shorelines.html

Kıyı şeridini düzenlemek için SVG'ye dönüştürdüm, daha sonra Inkscape ile düzenledim, sonra da lat / lon metin dosyasına ("MATLAB" formatına) dönüştürdüm.

Tüm Python kodları aşağıda yer almaktadır.

# ---------------------------------------------------------------
def plot_lines(mymap, lons, lats, **kwargs) :
    """Plots a custom coastline.  This plots simple lines, not
    ArcInfo-style SHAPE files.

    Args:
        lons: Longitude coordinates for line segments (degrees E)
        lats: Latitude coordinates for line segments (degrees N)

    Type Info:
        len(lons) == len(lats)
        A NaN in lons and lats signifies a new line segment.

    See:
        giss.noaa.drawcoastline_file()
    """

    # Project onto the map
    x, y = mymap(lons, lats)

    # BUG workaround: Basemap projects our NaN's to 1e30.
    x[x==1e30] = np.nan
    y[y==1e30] = np.nan

    # Plot projected line segments.
    mymap.plot(x, y, **kwargs)


# Read "Matlab" format files from NOAA Coastline Extractor.
# See: http://www.ngdc.noaa.gov/mgg/coast/

lineRE=re.compile('(.*?)\s+(.*)')
def read_coastline(fname, take_every=1) :
    nlines = 0
    xdata = array.array('d')
    ydata = array.array('d')
    for line in file(fname) :
#        if (nlines % 10000 == 0) :
#            print 'nlines = %d' % (nlines,)
        if (nlines % take_every == 0 or line[0:3] == 'nan') :
            match = lineRE.match(line)
            lon = float(match.group(1))
            lat = float(match.group(2))

            xdata.append(lon)
            ydata.append(lat)
        nlines = nlines + 1


    return (np.array(xdata),np.array(ydata))

def drawcoastline_file(mymap, fname, **kwargs) :
    """Reads and plots a coastline file.
    See:
        giss.basemap.drawcoastline()
        giss.basemap.read_coastline()
    """

    lons, lats = read_coastline(fname, take_every=1)
    return drawcoastline(mymap, lons, lats, **kwargs)
# =========================================================
# coastline2svg.py
#
import giss.io.noaa
import os
import numpy as np
import sys

svg_header = """<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->

<svg
   xmlns:dc="http://purl.org/dc/elements/1.1/"
   xmlns:cc="http://creativecommons.org/ns#"
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg"
   version="1.1"
   width="360"
   height="180"
   id="svg2">
  <defs
     id="defs4" />
  <metadata
     id="metadata7">
    <rdf:RDF>
      <cc:Work
         rdf:about="">
        <dc:format>image/svg+xml</dc:format>
        <dc:type
           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
        <dc:title></dc:title>
      </cc:Work>
    </rdf:RDF>
  </metadata>
  <g
     id="layer1">
"""

path_tpl = """
    <path
       d="%PATH%"
       id="%PATH_ID%"
       style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
"""

svg_footer = "</g></svg>"




# Set up paths
data_root = os.path.join(os.environ['HOME'], 'data')
#modelerc = giss.modele.read_modelerc()
#cmrun = modelerc['CMRUNDIR']
#savedisk = modelerc['SAVEDISK']

ifname = sys.argv[1]
ofname = ifname.replace('.dat', '.svg')

lons, lats = giss.io.noaa.read_coastline(ifname, 1)

out = open(ofname, 'w')
out.write(svg_header)

path_id = 1
points = []
for lon, lat in zip(lons, lats) :
    if np.isnan(lon) or np.isnan(lat) :
        # Process what we have
        if len(points) > 2 :
            out.write('\n<path d="')
            out.write('m %f,%f L' % (points[0][0], points[0][1]))
            for pt in points[1:] :
                out.write(' %f,%f' % pt)
            out.write('"\n   id="path%d"\n' % (path_id))
#            out.write('style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"')
            out.write(' />\n')
            path_id += 1
        points = []
    else :
        lon += 180
        lat = 180 - (lat + 90)
        points.append((lon, lat))


out.write(svg_footer)
out.close()

# =============================================================
# svg2coastline.py

import os
import sys
import re

# Reads the output of Inkscape's "Plain SVG" format, outputs in NOAA MATLAB coastline format

mainRE = re.compile(r'\s*d=".*"')
lineRE = re.compile(r'\s*d="(m|M)\s*(.*?)"')

fname = sys.argv[1]


lons = []
lats = []
for line in open(fname, 'r') :
    # Weed out extraneous lines in the SVG file
    match = mainRE.match(line)
    if match is None :
        continue

    match = lineRE.match(line)

    # Stop if something is wrong
    if match is None :
        sys.stderr.write(line)
        sys.exit(-1)

    type = match.group(1)[0]
    spairs = match.group(2).split(' ')
    x = 0
    y = 0
    for spair in spairs :
        if spair == 'L' :
            type = 'M'
            continue

        (sdelx, sdely) = spair.split(',')
        delx = float(sdelx)
        dely = float(sdely)
        if type == 'm' :
            x += delx
            y += dely
        else :
            x = delx
            y = dely
        lon = x - 180
        lat = 90 - y
        print '%f\t%f' % (lon, lat)
    print 'nan\tnan'

2
2018-01-30 20:48





Tamam sanırım kısmi bir çözümüm var.

Temel fikir, drawcoastlines () tarafından kullanılan yolların boyut / alan tarafından sıralanmasıdır. Demek ki ilk N yolu (çoğu uygulama için) ana kara kütleleri ve gölleri ve daha küçük adalar ve nehirler.

Sorun, istediğiniz ilk N yolunun projeksiyona (örneğin, küresel, kutupsal, bölgesel), eğer area_thresh uygulanmışsa ve göller ya da küçük adalar vb. İsteyip istemediğinize bağlı olacaktır. Başka bir deyişle, çimdiklemeniz gerekecektir. Bu uygulama başına.

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

mp = 'cyl'
m = Basemap(resolution='c',projection=mp,lon_0=0,area_thresh=200000)

fill_color = '0.9'

# If you don't want lakes set lake_color to fill_color
m.fillcontinents(color=fill_color,lake_color='white')

# Draw the coastlines, with a thin line and same color as the continent fill.
coasts = m.drawcoastlines(zorder=100,color=fill_color,linewidth=0.5)

# Exact the paths from coasts
coasts_paths = coasts.get_paths()

# In order to see which paths you want to retain or discard you'll need to plot them one
# at a time noting those that you want etc. 
for ipoly in xrange(len(coasts_paths)):
    print ipoly
    r = coasts_paths[ipoly]
    # Convert into lon/lat vertices
    polygon_vertices = [(vertex[0],vertex[1]) for (vertex,code) in
                        r.iter_segments(simplify=False)]
    px = [polygon_vertices[i][0] for i in xrange(len(polygon_vertices))]
    py = [polygon_vertices[i][1] for i in xrange(len(polygon_vertices))]
    m.plot(px,py,'k-',linewidth=1)
    plt.show()

Çizim yapmayı durdurmak için ilgili ipoliyi (poly_stop) öğrendikten sonra, böyle bir şey yapabilirsiniz ...

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

mproj = ['nplaea','cyl']
mp = mproj[0]

if mp == 'nplaea':
    m = Basemap(resolution='c',projection=mp,lon_0=0,boundinglat=30,area_thresh=200000,round=1)
    poly_stop = 10
else:
    m = Basemap(resolution='c',projection=mp,lon_0=0,area_thresh=200000)
    poly_stop = 18
fill_color = '0.9'

# If you don't want lakes set lake_color to fill_color
m.fillcontinents(color=fill_color,lake_color='white')

# Draw the coastlines, with a thin line and same color as the continent fill.
coasts = m.drawcoastlines(zorder=100,color=fill_color,linewidth=0.5)

# Exact the paths from coasts
coasts_paths = coasts.get_paths()

# In order to see which paths you want to retain or discard you'll need to plot them one
# at a time noting those that you want etc. 
for ipoly in xrange(len(coasts_paths)):
    if ipoly > poly_stop: continue
    r = coasts_paths[ipoly]
    # Convert into lon/lat vertices
    polygon_vertices = [(vertex[0],vertex[1]) for (vertex,code) in
                        r.iter_segments(simplify=False)]
    px = [polygon_vertices[i][0] for i in xrange(len(polygon_vertices))]
    py = [polygon_vertices[i][1] for i in xrange(len(polygon_vertices))]
    m.plot(px,py,'k-',linewidth=1)
plt.show()

enter image description here


1
2018-01-14 19:13





@ Sampo-smolander hakkındaki yorumuma göre

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

fig = plt.figure(figsize=(8, 4.5))
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.00)
m = Basemap(resolution='c',projection='robin',lon_0=0)
m.fillcontinents(color='gray',lake_color='white',zorder=2)
coasts = m.drawcoastlines(zorder=1,color='white',linewidth=0)
coasts_paths = coasts.get_paths()

ipolygons = range(83) + [84]
for ipoly in xrange(len(coasts_paths)):
    r = coasts_paths[ipoly]
    # Convert into lon/lat vertices
    polygon_vertices = [(vertex[0],vertex[1]) for (vertex,code) in
                        r.iter_segments(simplify=False)]
    px = [polygon_vertices[i][0] for i in xrange(len(polygon_vertices))]
    py = [polygon_vertices[i][1] for i in xrange(len(polygon_vertices))]
    if ipoly in ipolygons:
        m.plot(px,py,linewidth=0.5,zorder=3,color='black')
    else:
        m.plot(px,py,linewidth=0.5,zorder=4,color='grey')
plt.savefig('world2.png',dpi=100)

enter image description here


1
2018-01-24 17:34