Python: Recorta el cuadrado de la trama. El cuadrado tiene una longitud de lado de 30 kilómetros (18 millas) y tiene el centro en las coordenadas dadas

Question

El siguiente código muestra las coordenadas de una ciudad de Polonia indicadas por el usuario, por ejemplo Kock .

import subprocess
import sys
from pyproj import Proj, transform
import warnings

warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=DeprecationWarning)

name = input("Enter the name of the object (the town, e.g. Lublin or Kock): ") #Lublin, Kock

query = "SELECT ST_Centroid(ogr_geocode('{}, województwo lubelskie'))".format(name)
command = ["ogrinfo", ":memory:", "-q", "-sql", query]

proc = subprocess.run(command,
                      stdout=subprocess.PIPE,
                      encoding="utf-8")

if proc.returncode:
    sys.exit("Failed to execute ogrinfo command.")

lines = proc.stdout.split('\n')

found = False
for line in lines:
    line = line.strip()
    if line.startswith("POINT"):
        found = True

        _, x, y = line.split()
        x = x.lstrip("(")
        y = y.rstrip(")")
        print("")
        print("{}:\n".format(name))
        print("WGS84 (EPSG: 4326) coordinates")
        print("longitude:  ", x)
        print("latitude:  ", y)
        print("")
        break
if not found:
    sys.exit("City not found.")

inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:2180')
x1,y1 = x,y
x2,y2 = transform(inProj,outProj,x1,y1)
print("PUWG 1992 (EPSG: 2180) coordinates")
print ("X: ",x2)
print ("Y: ",y2)

Salida

Enter the name of the object (the town, e.g. Lublin or Kock): Kock

Kock:

WGS84 (EPSG: 4326) coordinates
longitude:   22.4458062648351
latitude:   51.654282346758

PUWG 1992 (EPSG: 2180) coordinates
X:  738269.852848907
Y:  426492.45745886303

En Salida, el script muestra las coordenadas de la ciudad en EPSG:4326, y después en EPSG: 2180.

Tengo un DEM para toda Polonia en el Polonia.tif en el sistema de coordenadas EPSG: 2180.

Cómo ampliar el script existente para que desde el Polonia.tif raster cortó un fragmento cuadrado (con una longitud de lado de 30 kilómetros [18 millas]) con el centro donde las coordenadas EPSG:2180?

Answer 1

Podrías llamar a gdal_translate en la trama de entrada y utilice la función -projwin para especificar los límites deseados.

Ya que conoce la posición de la ciudad x & y para calcular los lados de tu caja es simple matemática.

side = 30000 # 30km
h_side = side//2 # half a side (15km)
left = x - h_side
right = x + h_side
top = y + h_side
bottom = y - h_side

basta con llamar a gdal_translate -projwin left top right bottom ...