2

I have a code that will read a .tif and split the image into a specific number of rows and columns (code 1). I modified it to be able to be able to easily set the wanted image size (by pixels) (code 2). The wanted size wont be perfect because I am using the celi function to round to the nearest int. That is not an issue. Basically I am modifying QGIS way of creating a grid and clipping to each cell.

With the split images, I am running various classification methods to detect specific things. Unfortunately I cannot share my data with you. The problem that I am facing is that the left side of my images are not being split and therefor no images are able to run with the classification. You can see how in image 1 I have my classification areas and in image 2 the classification areas do not reach all the way over. There is a specific section where they stop. And yes, there are areas here that would be identified.

I've tried this method on a few images now and the same thing occurs. It's always on the left side. I have set my rows/columns to 3 and 3 and the missing section is no longer missing. It does split. So I do not understand why when I increase my columns and rows to a image size of say 600x600 (pixels), I end up missing data.

Is there something in my code that I am missing or does anyone have any suggestions?

# Setting the directory
os.chdir(r"D:\ortho")

Loading in the image


rimg = gdal.Open("image.tif")


Upper left corner of the minX and maxY


gt = rimg.GetGeoTransform()


xmin = gt[0]
ymax = gt[3]


res = gt[1]
xlen = res * rimg.RasterXSize # units is important UTM or WGS
ylen = res * rimg.RasterYSize


how many tiles you want to have in each row and column


xdiv = 70
ydiv = 105


Determining the size of each new image


xsize = xlen/xdiv 
ysize = ylen/ydiv


print(xsize)
print(ysize)


xsteps = [xmin + xsize * i for i in range(xdiv+1)] # plut because we start in the left top corner where X is at its lowest
ysteps = [ymax - ysize * i for i in range(ydiv+1)] # minus because we start in the left top corner where Y is at its highest


for i in range(xdiv):
    for j in range(ydiv):
        xmin = xsteps[i]
        xmax = xsteps[i+1]
        ymax = ysteps[j]
        ymin = ysteps[j+1]


    # Splices the image up into the set divs and saves them.
    gdal.Warp("D:/data/image" + str(i)+str(j) + ".tif",rimg,
              outputBounds = (xmin,ymin,xmax,ymax),dstNodata = -9999)



wanted_x = 2000
wanted_y = 2000
y_pixel_len = wanted_y * res
x_pixel_len = wanted_x * res


# how many tiles you want to have in each row and column
xdiv = math.ceil(xlen/x_pixel_len)
ydiv = math.ceil(ylen/y_pixel_len)


# Determining the size of each new image
xsize = math.ceil(xlen/xdiv)
ysize = math.ceil(ylen/ydiv)

image 1 enter image description here

image 2 enter image description here

Edits There are some warnings that occur when I run the code. I think these are just images that contain no data. So I end up getting a black image. I do have code that will delete those images. Not sure if this helps???

Warning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range.

More than 1000 errors or warnings have been reported. No more will be reported from now.

Binx
  • 1,350
  • 1
  • 9
  • 35
  • 1
    If you will consider a rasterio solution, have a look at this answer – user2856 Feb 13 '21 at 05:28
  • Take a look to this answer: https://gis.stackexchange.com/questions/59444/how-to-prevent-qgis-clipper-to-clamp-nodata-values – Gakovî Feb 13 '21 at 10:05
  • @user2856 I'll take a a look at this tomorrow. But it doesn't necessarily solve why my code above is causing an error like that. – Binx Feb 15 '21 at 02:49
  • @XiJin looks like that is set within QGIS, how do I set it in my python code? – Binx Feb 15 '21 at 02:50
  • Yeah, that's why I didn't add as an answer. – user2856 Feb 15 '21 at 04:21
  • Try explicitly setting the target resolution. But without having the data it is just guesswork. Also have you zoomed in to the part with missing data? Sometimes interpolation and zoom level makes it seem like nothing is there. – JonasV Feb 16 '21 at 11:04
  • @JonasV I have been working on trying to find which resolution starts cutting out the images. Unfortunately the process of classifying the image takes a while, so the process is slow. I have zoomed in and there is no data in that region that is missing. – Binx Feb 16 '21 at 18:25
  • @user2856 Still working on your suggestion. I will have to say, I think it splits the images faster which is always good. – Binx Feb 16 '21 at 18:26
  • What else is happening with the image? Potentially the classification is the culprit – JonasV Feb 18 '21 at 09:42
  • @JonasV I am using the Dzetsaka: classification tool plugin to identify a vegetation species. I've sifted through the split images and cannot find the section that is missing. In order to classify this section I have to clip it out first and then split it. Just another step that should be avoidable. The clipped out section runs fine with the classification, results are produced. – Binx Feb 18 '21 at 16:59
  • @user2856 I was able to run your code through your link you provided and it did work. Unfortunately (should be asked in a different question), that code cannot be used in a batch file. A DLL failed to load model error on from rasterio._base import gdal_version is the culprit. I've seen quite a few github and stack questions regarding this issue too. – Binx Feb 19 '21 at 01:23
  • Anyways, your code works, but in my case I cannot use it because of the DLL failed to load error. Ultimately I need to automate the splitting and classification process. The code I posted above does work with the batch file, but as we can see is missing data. – Binx Feb 19 '21 at 01:25

0 Answers0