Python script, cutting area with gdal_translate using Spyder

1217
1
06-07-2020 07:07 AM
RafaelMichael
Occasional Contributor

I am attaching a file lubelskie.tiff (DEM in the reference system EPSG:2180, for one of the counties in Poland).

I need to write a python script, in which the user is asked to enter the name of the town belonging to that county (for example cities Lublin, Frampol, Parczew), and this script then cuts out from the lubelskie.tiff file a square with a side of 30 kilometers, whose center is in the center of the entered town. The cut area is to be saved to the city.tiff file, where the word "city" is to be replaced by the actual city name.

You can use the ogrinfo command (from the attached file zad_07.py) to find the coordinates of the city (in the WGS84 coordinate system (longitude, latitude)), and the gdaltransform command (from the attached file zad_08.py) to convert the WGS84 coordinates to coordinates EPSG: 2180. In addition, you must add gdal_translate command to cut the area.

The gdal_translate command can probably be used in the following way:

gdal_translate -projwin ulx uly lrx lry lubelskie.tif city.tif

where ulx, uly, lrx, lry are respectively: the x coordinate of the upper left vertex, y coordinate of the upper left vertex, the x coordinate of the bottom right vertex, the y coordinate of the lower right vertex of the cut out fragment.

The files are packed here: lubelskie.tiff, zad_7.py, zad_8.py.

https://www97.zippyshare.com/v/zyQXpmmx/file.html

Can you help me write this script? I am stuck and need help

0 Kudos
1 Reply
DavidPike
MVP Frequent Contributor

What part do you have a problem with (I can't help with any GDAL stuff unfortunately).

I would however:

1.Create a python list or dictionary of your Cities as a Key, and the WGS84 Lat Long as a value.

2. Run these through GDAL translate to get the projected coordinates

3. Do some simple maths to get your UL/BR coordinates for each city 30km square.

MyDictionary = {"Lublin" : [290004 370001, 320004 34001], "Frampol" : [ ...

workspace = r"C:\\.."

Input = "Lublin"

#create output path

outpath = Input+".tif"

outpath = os.path.join(workspace, outpath)

GDAL_function(MyDictionary[Input], outpath)

0 Kudos