Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
q_350771_68792.ipynb
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@hibajamal

This comment has been minimized.

Copy link

hibajamal commented May 28, 2020

What could be the reason for getting this output when the input image has been changed, and the bands kept the same:
image

@jdbcode

This comment has been minimized.

Copy link
Owner Author

jdbcode commented May 28, 2020

Did you also update the aoi geometry? It might be that the region of pixels that you are grabbing (based on the aoi) are outside of the image or its valid pixels. (?)

@hibajamal

This comment has been minimized.

Copy link

hibajamal commented May 28, 2020

Thanks so much for replying. Yeah, I did. I updated it for a specific lat and long value, and my region is a square on that lat, long (attaching code). Does it have something to do with the bands? Sorry for the questions rookie here!
image

@jdbcode

This comment has been minimized.

Copy link
Owner Author

jdbcode commented May 28, 2020

Ahhh, I think I see the problem. My example uses the SR image product and you are using TOA. The SR images are scaled by 10,000, so you need to change the line that defines rgb_img_test to:

rgb_img_test = (255*(rgb_img/0.35)).astype('uint8')
@hibajamal

This comment has been minimized.

Copy link

hibajamal commented May 28, 2020

Got it, it worked. Thank you so much!!!

@jdbcode

This comment has been minimized.

Copy link
Owner Author

jdbcode commented May 28, 2020

Happy to hear!

@hibajamal

This comment has been minimized.

Copy link

hibajamal commented Jun 4, 2020

Hey, incredibly sorry for bothering you again but, is there any way to set a scale or maxPixels parameter during (or after) the extraction of band values in sampleRectangle()?
It's just that I for a couple of regions and images, I keep getting the too many pixels error, and am out of ways to fix it. Help please (again, sorry).
image

@jdbcode

This comment has been minimized.

Copy link
Owner Author

jdbcode commented Jun 4, 2020

Hi @hibajamal

That is a limitation of the function. The only way to get around it is to make a smaller rectangle or resample the image to have larger pixels like this:

band_arrs = img.reproject(crs=img.projection().crs(), scale=100).sampleRectangle(region=aoi)
@hibajamal

This comment has been minimized.

Copy link

hibajamal commented Jun 4, 2020

Thanks so much for the response, I tried limiting the size by fixing a small rectangle and now this, adding a high scale but that leads to giving a completely white image, or if the rectangle is too small, then a black image (since i kept defaultValue=0 in sampleRectangle).

@hibajamal

This comment has been minimized.

Copy link

hibajamal commented Jun 5, 2020

Hey, @jdbcode,
Also, will rgb_img_test have to be changed or scaled according to the scale in img.reproject or not? Sorry I am just trying figure out why the output image looks like this (it is also the median of an image collection if that matters):
image

@oscmansan

This comment has been minimized.

Copy link

oscmansan commented Jun 25, 2020

Thanks for the gist! Do you have any idea why if I use the COPERNICUS/S2_SR dataset sampleRectangle() returns a single pixel? I am sampling a region of 1000x1000 m and the resolution of the bands is 10 m, so I was expecting 100x100 pixels.

Here's my code:

collection = ee.ImageCollection('COPERNICUS/S2_SR')
collection = collection.filterDate('2020-01-01', '2020-01-30')
collection = collection.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
collection = collection.map(maskS2clouds)

image = collection.mean()

region = ee.Geometry.Point(83.277, 17.7009).buffer(1000).bounds()
patch = image.sampleRectangle(region)
plt.imshow(np.dstack(([patch.get(b).getInfo() for b in ['B4', 'B3', 'B2']])))
@jdbcode

This comment has been minimized.

Copy link
Owner Author

jdbcode commented Jun 26, 2020

Image composites (collection.mean()) have a default projection of WGS84 with 1-degree scale and since the region is less than 1 degree, only a single pixel is returned. Set the projection of the mean composite manually with setDefaultProjection(). See suggested code edits for possible solution (untested).

region = ee.Geometry.Point(83.277, 17.7009).buffer(1000).bounds() # <-- move above collection def
collection = ee.ImageCollection('COPERNICUS/S2_SR')
collection = collection.filterBounds(region) # <-- filter by region intersection
collection = collection.filterDate('2020-01-01', '2020-01-30')
collection = collection.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
collection = collection.map(maskS2clouds)

proj = collection.first().projection() # <-- get the projection of an image
image = collection.mean().setDefaultProjection(proj); # <--set projection of the composite

patch = image.sampleRectangle(region)
plt.imshow(np.dstack(([patch.get(b).getInfo() for b in ['B4', 'B3', 'B2']])))
@oscmansan

This comment has been minimized.

Copy link

oscmansan commented Jun 27, 2020

Thanks a lot @jdbcode! I didn't know about the default projection of image composites. Actually, I don't even need a composite, I just want to sample an image from that region (e.g. the most recent pass). So I fixed my code as follows:

collection = ee.ImageCollection('COPERNICUS/S2_SR')
collection = collection.filterDate('2017-03-28', '2020-06-27')
collection = collection.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
collection = collection.map(maskS2clouds)

region = ee.Geometry.Point(7.411762, 43.725962).buffer(1000).bounds()
image = collection.filterBounds(region).limit(1, 'system:index', False).first()
bands = image.sampleRectangle(region, defaultValue=0)
patch = np.dstack(([bands.get(b).getInfo() for b in ['B4', 'B3', 'B2']]))

plt.imshow(patch/0.3)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.