0

I have a large sentinel-2 image and I want to clip it according polygons which are stored in a geojson file.

My code is:

df = gpd.read_file(os.path.join(geojson_folder,geojson_file))
df.head()

    id  name    geometry
0   1   almond  (POLYGON ((749159.7552032885 3643270.487528285...
1   2   olive   (POLYGON ((749639.4218857342 3643450.362534203...
2   3   banana  (POLYGON ((749055.617041968 3643819.579651611,...
3   4   orange  (POLYGON ((748500.2135149258 3643882.693688774...
4   5   apple   (POLYGON ((747988.9898138981 3643119.013839092...

Where gpd is geopandas, and the crs is espg:32636

To perform clip (using rasterio.mask) for each polygon by each sentinel-2 raster (format jp2) I use:

os.chdir(output)

coords = [json.loads(df.to_json())['features'][i]['geometry'] for i in range((df.shape[0]))]
field_name = [crop for crop in df.name.tolist()]
for poly,fn in zip(coords,field_name):
    os.mkdir(fn)
    for sn2_folder in os.listdir(src):
        if sn2_folder.endswith('.zip'):
            pass
        else:
            for sn2_band in os.listdir(os.path.join(src,sn2_folder)):

                data = rasterio.open(os.path.join(src,sn2_folder,sn2_band))
                epsg_code = int(data.crs.data['init'][5:])
                srs = osr.SpatialReference()
                srs.ImportFromEPSG(epsg_code)

                try:            
                    out_img, out_transform = mask(dataset=data, shapes=[poly], crop=True, all_touched=False)
                    out_meta = data.meta.copy()
                    new_name = sn2_band
                    out_tif = os.path.join(fn,new_name)
                    out_meta.update({"driver": "JP2OpenJPEG",
                         "height": out_img.shape[1],
                         "width": out_img.shape[2],
                         "transform": out_transform,
                         "crs": srs.ExportToProj4()}
                                 )

                    with rasterio.open(out_tif, "w", **out_meta) as dest:
                        dest.write(out_img)
                except:
                    pass

print ("finish")

The outcome seems satisfying. I get clipped rasters, which are correctly geo-referenced, with the desired pixels size etc. However, I've noticed that the pixel values of the original raster are different from the clipped raster. Why is that happening? is there some hidden resample (I didn't see anything in the official documents)?

UPDATE: After playing with this some more, the (weird) solution is to change the "JP2OpenJPEG" to "Gtiff". For some reason this creates jp2 files with matching values between original and clipped rasters.

  • 1
    How different are you talking? JPEG 2000 is a lossy format so it'll never be able to represent the original values exactly – mikewatt Mar 19 at 19:56
  • 1
    JPEG 2000 has also a lossless mode but probably it is not used for those images. – user30184 Mar 19 at 21:35
  • @mikewatt- let's say the original value is 902, then the new value is 870, differences in that magnitude. The differences are not huge, I was just wondering why the difference exists in the first place, it is a simple process, maybe is has something to do with geopandas? – user88484 Mar 20 at 4:22
  • @user88484 highly unlikely that it has anything to do with geopandas. You can tests this yourself by using an array mask of your own creation – Paul H Mar 26 at 22:42

Your Answer

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

Browse other questions tagged or ask your own question.