3

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
Improve this question
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 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/… Commented Mar 8, 2024 at 9:43

2 Answers 2

Reset to default
1

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.

Improve this answer
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. Commented Jul 18, 2022 at 11:54
-1

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

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.