2

I am trying to clip a raster image into a small image using coordinates in a CSV file. In the CSV file, I have minx,maxx,miny,maxy as the coordinates which makes the bounding box and clips the images.

I have a total of 100,000 coordinates and it will create 50,000 images. I am trying to number the images in ascending numbers. Example 0,1,2,3.... till 50,000. The problem is that when clipping it's escaping some images and only clipping 26,000 images. I just cannot find the bug in my code.

I thought it was the coordinate but when checking the coordinates it seems fine, the bug somehow is in my code.

import pandas as pd 
from osgeo import gdal

df = pd.read_csv("G:\\satellite\\align\\satellite.txt")
fs = gdal.Open('G:\\satellite\\align\\satellite.tif')

i, j = 0, 0

while True:
    minx = df.iloc[i][0]
    miny = df.iloc[i][1]
    maxx = df.iloc[i+1][0]
    maxy = df.iloc[i+1][1]
    
    out_tif = r"G:\\satellite\\training_images\\satellite_rename\\{0}.tif".format(j)
    ds = gdal.Translate(out_tif, fs, projWin = [minx, miny, maxx, maxy])
    
    j = j+1
    i += 2
    if i > 100000: # The iteration depends on how many co-ordinates you have
        break
4
  • 2
    are your coordinates and rasters in the same projection? Commented Jul 12, 2022 at 21:39
  • Yes i took the coordinates from point feature class and they are in the same projection
    – user123
    Commented Jul 13, 2022 at 7:43
  • another question. have you stopped the process after one clip to view it? along with the coordinates. does the clip make sense? Commented Jul 13, 2022 at 19:10
  • why not host the Geotiff file in a geoserver as a wms service ? That way you can access parts of it using any desired bounding box. Check this: stackoverflow.com/questions/3015079/…
    – Abhilshit
    Commented Mar 8 at 9:43

2 Answers 2

0

According to gdal python docs:

projWin --- subwindow in projected coordinates to extract: [ulx, uly, lrx, lry]

uly would be 'upper-left' which corresponds to a maximum for y. Try the following for your projWin:

projWin = [minx, maxy, maxx, miny]

It could also explain why your script stops making crops. If the raster is not square (they rarely are) then at some point in your iteration the coordinates may be out of bounds, so to speak, and thus no crop.

1
  • I tired this and it directly throws error. One weird thing is i also have other images and the same numbers are not being cropped.
    – user123
    Commented Jul 18, 2022 at 11:54
0

You can clip using rasterio.mask and geopandas:

import geopandas as gpd
import rasterio, os
from shapely import box
from rasterio.mask import mask

raster_file = r"/home/bera/Drives/data1/gisdata/Sentinel2/Relativt_molnfria_spridda_bilder/Lov_av/2020/S2B_MSIL2A_20200327T101629_N0214_R065_T33VVD_20200327T134849.SAFE/GRANULE/L2A_T33VVD_A015964_20200327T101842/IMG_DATA/R10m/T33VVD_20200327T101629_B08_10m.jp2"
df = gpd.pd.read_csv(r"/home/bera/Desktop/gistest/satellite/coordinates.csv")
out_folder = r"/home/bera/Desktop/gistest/satellite/clipped_images/"

# df.head(1)
#    id    minx    maxx     miny     maxy
# 0   1  399960  402960  6396020  6400020

#Create a geodataframe with polygon geometries
df["geometry"] = df.apply(lambda x: box(xmin=x.minx, ymin=x.miny, xmax=x.maxx, ymax=x.maxy), axis=1)
df = gpd.GeoDataFrame(data=df, geometry=df.geometry, crs=32633)

#Clip the raster with each polygon
source = rasterio.open(raster_file)
shapes = df.geometry.tolist()

for num, shape in enumerate(shapes, 1):
    out_file = os.path.join(out_folder, f"clip_{num}.tif") #Name the output clip_1.tif, clip_2.tif, ...
    out_image, out_transform = mask(dataset=source, shapes=[shape], crop=True) #Clip the raster
    with rasterio.open(out_file, mode="w", driver="GTiff", width=out_image.shape[2], #Write to file
                       height=out_image.shape[1], count=1, crs=source.crs, transform=out_transform,
                       dtype=out_image.dtype, compression="lzw") as destination:
        destination.write(out_image)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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