1

I have a GRID-file with a quite large extent. For an automatic workflow I need the maximum extent/bounding-box for all z-values that are over a certain range. For example: Get the maximum extent for all z-values that are greater or equal than 0.05

Is there a possibility to do that with Python, GDAL or GMT-tools (generic mapping tools), so that it can run automatically?

1

not the best solution but solved the problem for me. For anybody who wants to know how I did it:

gdal_contour -i 1000 -off 0.05 input.grd output.shp

-i = some very high number -off = minimum pixel value

then:

ogrinfo -al -so output.shp

and parse the ogrinfo output to get the extent coords

EDIT Better Solution (Python3 with Subprocess):

gdal_calc.py -A input.grd --oufile=temp.tif --calc=logical_an(A>=min, A<=max)
gdal_polygonize.py temp.tif -f GMT temp.gmt

with Python3:

gmt_file = open(temp.gmt)
extent_temp = gmt_file.readlines[1]
extent = re.findall("(-?\d+.\d+)",extent_temp)

west = extent_w_time[0]   
east = extent_w_time[1]   
south = extent_w_time[2]   
north = extent_w_time[3]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.