Profil wysokości 10 km z każdej strony linii

Jak mogę uzyskać profil wysokości dla pasma terenu?

Należy wziąć pod uwagę najwyższe wzniesienie w promieniu 10 km (po każdej stronie zdefiniowanej linii).

Mam nadzieję, że moje pytanie jest jasne. Z góry bardzo dziękuję.

Kontynuując komentarze, oto wersja, która działa z prostopadłymi segmentami linii. Proszę używać ostrożnie, ponieważ nie przetestowałem go dokładnie!

Ta metoda jest o wiele bardziej niezręczna niż odpowiedź @ whubera - częściowo dlatego, że nie jestem bardzo dobrym programistą, a częściowo dlatego, że przetwarzanie wektorowe jest trochę kłopotliwe. Mam nadzieję, że przynajmniej zaczniesz, jeśli potrzebne są prostopadłe segmenty linii.

Musisz uruchomić pakiety Shapely , Fiona i Numpy Python (wraz z ich zależnościami), aby to uruchomić.

#-------------------------------------------------------------------------------
# Name:        perp_lines.py
# Purpose:     Generates multiple profile lines perpendicular to an input line
#
# Author:      JamesS
#
# Created:     13/02/2013
#-------------------------------------------------------------------------------
""" Takes a shapefile containing a single line as input. Generates lines
    perpendicular to the original with the specified length and spacing and
    writes them to a new shapefile.

    The data should be in a projected co-ordinate system.
"""

import numpy as np
from fiona import collection
from shapely.geometry import LineString, MultiLineString

# ##############################################################################
# User input

# Input shapefile. Must be a single, simple line, in projected co-ordinates
in_shp = r'D:\Perp_Lines\Centre_Line.shp'

# The shapefile to which the perpendicular lines will be written
out_shp = r'D:\Perp_Lines\Output.shp'

# Profile spacing. The distance at which to space the perpendicular profiles
# In the same units as the original shapefile (e.g. metres)
spc = 100

# Length of cross-sections to calculate either side of central line
# i.e. the total length will be twice the value entered here.
# In the same co-ordinates as the original shapefile
sect_len = 1000
# ##############################################################################

# Open the shapefile and get the data
source = collection(in_shp, "r")
data = source.next()['geometry']
line = LineString(data['coordinates'])

# Define a schema for the output features. Add a new field called 'Dist'
# to uniquely identify each profile
schema = source.schema.copy()
schema['properties']['Dist'] = 'float'

# Open a new sink for the output features, using the same format driver
# and coordinate reference system as the source.
sink = collection(out_shp, "w", driver=source.driver, schema=schema,
                  crs=source.crs)

# Calculate the number of profiles to generate
n_prof = int(line.length/spc)

# Start iterating along the line
for prof in range(1, n_prof+1):
    # Get the start, mid and end points for this segment
    seg_st = line.interpolate((prof-1)*spc)
    seg_mid = line.interpolate((prof-0.5)*spc)
    seg_end = line.interpolate(prof*spc)

    # Get a displacement vector for this segment
    vec = np.array([[seg_end.x - seg_st.x,], [seg_end.y - seg_st.y,]])

    # Rotate the vector 90 deg clockwise and 90 deg counter clockwise
    rot_anti = np.array([[0, -1], [1, 0]])
    rot_clock = np.array([[0, 1], [-1, 0]])
    vec_anti = np.dot(rot_anti, vec)
    vec_clock = np.dot(rot_clock, vec)

    # Normalise the perpendicular vectors
    len_anti = ((vec_anti**2).sum())**0.5
    vec_anti = vec_anti/len_anti
    len_clock = ((vec_clock**2).sum())**0.5
    vec_clock = vec_clock/len_clock

    # Scale them up to the profile length
    vec_anti = vec_anti*sect_len
    vec_clock = vec_clock*sect_len

    # Calculate displacements from midpoint
    prof_st = (seg_mid.x + float(vec_anti[0]), seg_mid.y + float(vec_anti[1]))
    prof_end = (seg_mid.x + float(vec_clock[0]), seg_mid.y + float(vec_clock[1]))

    # Write to output
    rec = {'geometry':{'type':'LineString', 'coordinates':(prof_st, prof_end)},
           'properties':{'Id':0, 'Dist':(prof-0.5)*spc}}
    sink.write(rec)

# Tidy up
source.close()
sink.close()

Poniższy obraz pokazuje przykład danych wyjściowych skryptu. Wpisujesz plik kształtu reprezentujący linię środkową i określasz długość linii prostopadłych oraz ich odstępy. Dane wyjściowe to nowy plik kształtu zawierający czerwone linie na tym obrazie, z których każdy ma powiązany atrybut określający jego odległość od początku profilu.

Jak powiedział @whuber w komentarzach, po przejściu do tego etapu reszta jest dość łatwa. Poniższy obraz pokazuje inny przykład z wyjściem dodanym do ArcMap.

Użyj narzędzia Cecha do rastra, aby przekształcić linie prostopadłe w kategoryczny raster. Ustaw raster VALUEna Distpole wyjściowego pliku kształtów. Należy również pamiętać, aby ustawić narzędzie Environmentstak, że Extent, Cell sizei Snap rastersą takie same jak dla podstawowego DEM. Powinieneś otrzymać rastrową reprezentację swoich linii, mniej więcej tak:

Na koniec przekonwertuj ten raster na siatkę całkowitą (za pomocą narzędzia Int lub kalkulatora rastrowego) i użyj go jako stref wejściowych dla statystyki strefowej jako narzędzia tabeli . Powinieneś otrzymać tabelę wyników taką jak ta:

VALUEPola w tabeli daje dystans od początku oryginalnej linii profilu. Pozostałe kolumny podają różne statystyki (maksimum, średnią itp.) Dla wartości w każdym transecie. Możesz użyć tej tabeli do wykreślenia swojego profilu podsumowującego.

Uwaga: Jednym oczywistym problemem związanym z tą metodą jest to, że jeśli twoja pierwotna linia jest bardzo niepewna, niektóre linie przecięcia mogą się nakładać. Narzędzia statystyk strefowych w ArcGIS nie mogą poradzić sobie z nakładającymi się strefami, więc gdy tak się stanie, jedna z linii transektu będzie miała pierwszeństwo przed drugą. To może, ale nie musi stanowić problemu w tym, co robisz.

Powodzenia!

Najwyższe wzniesienie w promieniu 10 km to maksymalna wartość sąsiedztwa obliczona z okrągłym promieniem 10 km, więc po prostu wyodrębnij profil maksymalnej siatki sąsiedztwa wzdłuż trajektorii.

Przykład

Oto cieniowany wzgórzem DEM z trajektorią (czarna linia biegnąca od dołu do góry):

Ten obraz ma wymiary około 17 na 10 kilometrów. Aby zilustrować tę metodę, wybrałem promień zaledwie 1 km zamiast 10 km. Bufor o długości 1 km jest zaznaczony na żółto.

Maksymalne sąsiedztwo DEM zawsze będzie wyglądać trochę dziwnie, ponieważ będzie miało tendencję do zwiększania wartości w punktach, w których jedno maksimum (być może szczyt) spada nieco ponad 10 km, a drugie maksimum na innej wysokości znajduje się w odległości 10 km . W szczególności szczyty wzgórz, które dominują w ich otoczeniu, przyczynią się do powstania idealnych kręgów wartości wyśrodkowanych w punkcie maksymalnej wysokości lokalnej:

Ciemniej jest wyżej na tej mapie.

Oto wykres profili oryginalnego DEM (niebieski) i maksimum sąsiedztwa (czerwony):

Obliczono go, dzieląc trajektorię na regularnie rozmieszczone punkty w odległości 0,1 km (zaczynając od południowego końca), wyodrębniając wzniesienia w tych punktach i wykonując połączony wykres rozproszenia powstałych potrójnych (odległość od początku, wzniesienia, wzniesienia maksymalnego). Odległość między punktami wynosząca 0,1 km została wybrana, aby była znacznie mniejsza niż promień bufora, ale wystarczająco duża, aby obliczenia przebiegały szybko (było to natychmiastowe).

Miałem ten sam problem i wypróbowałem rozwiązanie Jamesa S., ale nie mogłem zmusić GDAL do współpracy z Fioną.

Następnie odkryłem algorytm SAGA „Profile krzyżowe” w QGIS 2.4 i uzyskałem dokładnie taki wynik, jaki chciałem i zakładam, że też go szukasz (patrz poniżej).

Dla każdego, kto jest zainteresowany, oto zmodyfikowana wersja kodu JamesS tworząca linie prostopadłe przy użyciu tylko bibliotek numpy i osgeo. Dzięki JamesS jego odpowiedź bardzo mi dziś pomogła!

import osgeo
from osgeo import ogr
import numpy as np

# ##############################################################################
# User input

# Input shapefile. Must be a single, simple line, in projected co-ordinates
in_shp = r'S:\line_utm_new.shp'

# The shapefile to which the perpendicular lines will be written
out_shp = r'S:\line_utm_neu_perp.shp'

# Profile spacing. The distance at which to space the perpendicular profiles
# In the same units as the original shapefile (e.g. metres)
spc = 100

# Length of cross-sections to calculate either side of central line
# i.e. the total length will be twice the value entered here.
# In the same co-ordinates as the original shapefile
sect_len = 1000
# ##############################################################################

# Open the shapefile and get the data
driverShp = ogr.GetDriverByName('ESRI Shapefile')
sourceShp = driverShp.Open(in_shp, 0)
layerIn = sourceShp.GetLayer()
layerRef = layerIn.GetSpatialRef()

# Go to first (and only) feature
layerIn.ResetReading()
featureIn = layerIn.GetNextFeature()
geomIn = featureIn.GetGeometryRef()

# Define a shp for the output features. Add a new field called 'M100' where the z-value 
# of the line is stored to uniquely identify each profile
outShp = driverShp.CreateDataSource(out_shp)
layerOut = outShp.CreateLayer('line_utm_neu_perp', layerRef, osgeo.ogr.wkbLineString)
layerDefn = layerOut.GetLayerDefn() # gets parameters of the current shapefile
layerOut.CreateField(ogr.FieldDefn('M100', ogr.OFTReal))

# Calculate the number of profiles/perpendicular lines to generate
n_prof = int(geomIn.Length()/spc)

# Define rotation vectors
rot_anti = np.array([[0, -1], [1, 0]])
rot_clock = np.array([[0, 1], [-1, 0]])

# Start iterating along the line
for prof in range(1, n_prof):
    # Get the start, mid and end points for this segment
    seg_st = geomIn.GetPoint(prof-1) # (x, y, z)
    seg_mid = geomIn.GetPoint(prof)
    seg_end = geomIn.GetPoint(prof+1)

    # Get a displacement vector for this segment
    vec = np.array([[seg_end[0] - seg_st[0],], [seg_end[1] - seg_st[1],]])    

    # Rotate the vector 90 deg clockwise and 90 deg counter clockwise
    vec_anti = np.dot(rot_anti, vec)
    vec_clock = np.dot(rot_clock, vec)

    # Normalise the perpendicular vectors
    len_anti = ((vec_anti**2).sum())**0.5
    vec_anti = vec_anti/len_anti
    len_clock = ((vec_clock**2).sum())**0.5
    vec_clock = vec_clock/len_clock

    # Scale them up to the profile length
    vec_anti = vec_anti*sect_len
    vec_clock = vec_clock*sect_len

    # Calculate displacements from midpoint
    prof_st = (seg_mid[0] + float(vec_anti[0]), seg_mid[1] + float(vec_anti[1]))
    prof_end = (seg_mid[0] + float(vec_clock[0]), seg_mid[1] + float(vec_clock[1]))

    # Write to output
    geomLine = ogr.Geometry(ogr.wkbLineString)
    geomLine.AddPoint(prof_st[0],prof_st[1])
    geomLine.AddPoint(prof_end[0],prof_end[1])
    featureLine = ogr.Feature(layerDefn)
    featureLine.SetGeometry(geomLine)
    featureLine.SetFID(prof)
    featureLine.SetField('M100',round(seg_mid[2],1))
    layerOut.CreateFeature(featureLine)

# Tidy up
outShp.Destroy()
sourceShp.Destroy()