GeoscienceAustralia / dea-notebooks

Repository for Digital Earth Australia Jupyter Notebooks: tools and workflows for geospatial analysis with Open Data Cube and Xarray
https://docs.dea.ga.gov.au/notebooks/
Apache License 2.0
439 stars 127 forks source link

xarray Time #1164

Closed UQ-Raja-Ram-Aryal closed 5 months ago

UQ-Raja-Ram-Aryal commented 8 months ago
# Create a reusable query
query = {
    'y': lat_range,
    'x': lon_range,
    'time': time_range,
    'measurements': ['nbart_red', 'nbart_green', 'nbart_blue', 'nbart_nir_1','nbart_swir_3',
                     'satellite_azimuth','solar_azimuth','solar_zenith','satellite_view','relative_slope'],
    'resolution': (-10, 10),
    'output_crs': 'epsg:3577',
    'group_by': 'solar_day'
}   # changed res from 20 to 10, added 'nbart_swir_3'

# Load available data from Sentinel
ds = load_ard(dc=dc,
              products=['ga_s2am_ard_3', 'ga_s2bm_ard_3'],
              min_gooddata=0.90,
              **query)
print(ds) 

When I am running this code, I am getting Dimensions: time: 332 y: 15 x: 15 I expect my dataset in the 5-day or 10-day interval, but I am getting my Time value within 2 days/3 days/ 5 days/ 10 days intervals, why? I am interested in the exact date of the satellite visit with the same time interval. How I can get that time value?

robbibt commented 7 months ago

Hey @UQ-Raja-Ram-Aryal, can you let me know what values lat_range, lon_range, time_range have?

UQ-Raja-Ram-Aryal commented 7 months ago

Define Site Name

site_name = "Boyagin_Wandoo"

Function to rotate a point about another point by an angle theta (in degrees)

def rotate_point(x, y, x_origin, y_origin, theta): theta_rad = math.radians(theta) x_rotated = math.cos(theta_rad) (x - x_origin) - math.sin(theta_rad) (y - y_origin) + x_origin y_rotated = math.sin(theta_rad) (x - x_origin) + math.cos(theta_rad) (y - y_origin) + y_origin return x_rotated, y_rotated

Your original coordinates and ranges

lat = -32.4771 ### Accurate latitude of the supersites lon = 116.9386 ### Accurae Longitude of the supersites buffer_distance = 120 # in meters extension_distance = 230 # in meters

Convert distances to degrees

buffer_distance_deg = buffer_distance / 111000 extension_distance_deg = extension_distance / 111000

Create original corners of the box

lat1, lon1 = lat - extension_distance_deg, lon + extension_distance_deg lat2, lon2 = lat1, lon1 + buffer_distance_deg lat3, lon3 = lat1 + buffer_distance_deg, lon1 lat4, lon4 = lat1 + buffer_distance_deg, lon1 + buffer_distance_deg

Rotate by an angle theta

theta = 115 # Degrees lat_origin, lon_origin = lat, lon # Rotation origin

lat1, lon1 = rotate_point(lat1, lon1, lat_origin, lon_origin, theta) lat2, lon2 = rotate_point(lat2, lon2, lat_origin, lon_origin, theta) lat3, lon3 = rotate_point(lat3, lon3, lat_origin, lon_origin, theta) lat4, lon4 = rotate_point(lat4, lon4, lat_origin, lon_origin, theta)

Calculate new lat and lon ranges

lat_range = (min(lat1, lat2, lat3, lat4), max(lat1, lat2, lat3, lat4)) lon_range = (min(lon1, lon2, lon3, lon4), max(lon1, lon2, lon3, lon4)) These are the input variables

robbibt commented 7 months ago

@UQ-Raja-Ram-Aryal Thankyou - can you also share the actual values you end up with, e.g.

print(lat_range) 
print(lon_range) 
print(time_range) 
UQ-Raja-Ram-Aryal commented 7 months ago

Latitude range in meters: 159.47112585420342 Longitude range in meters: 159.47112585341472 (-32.47953891636641, -32.47810223955691) (116.93538948484118, 116.93682616165067) ('2015-01-01', '2023-12-31')

robbibt commented 7 months ago

Hey @UQ-Raja-Ram-Aryal , looking at the result of your query, I think what's going on is that you're lucky enough to be getting more observations than usual because your site is in the overlap between two satellite overpasses. For example, I ran:

ds.time.to_dataframe().tail(10)

To get a list of some of the most recent overpasses: image

Then went to DEA Maps to look at the satellite imagery from the 22nd and 24th of December:

image image

You can see that you're getting data from the east part of the overpass on the 22nd, and the west side of the overpass on the 24th.

This will vary site-by-site, so at some locations you might only get the standard ~5-day revisits, and at other sites like this it might be a bit more frequent!

robbibt commented 7 months ago

Also, out of curiosity: I noticed that you're loading 'satellite_azimuth','solar_azimuth','solar_zenith','satellite_view','relative_slope'. I'd love to know what you're doing with those bands! We often don't know exactly what parts of our data are being used for what, so it's always useful to hear about possible new use cases. 🙂

UQ-Raja-Ram-Aryal commented 7 months ago

I am working on monitoring and evaluating the Land Surface Phenology of Australian vegetation using Sentinel-2. So, the Time is more valuable for this work. In the same location how it will visit with two days. Those dates are in AM, why? So, How Can I adjust my data with a 5-day interval?

UQ-Raja-Ram-Aryal commented 7 months ago

Also, out of curiosity: I noticed that you're loading 'satellite_azimuth','solar_azimuth','solar_zenith','satellite_view','relative_slope'. I'd love to know what you're doing with those bands! We often don't know exactly what parts of our data are being used for what, so it's always useful to hear about possible new use cases. 🙂

I am also thinking of looking at the impact of the solar view angle information to track the vegetation productivity of the site. Which I am not going to use anymore. For those objectives, I am going to use Radiative Transfer Model.

robbibt commented 7 months ago

I am working on monitoring and evaluating the Land Surface Phenology of Australian vegetation using Sentinel-2. So, the Time is more valuable for this work. In the same location how it will visit with two days. Those dates are in AM, why? So, How Can I adjust my data with a 5-day interval?

The times you get on data loaded from the datacube are in UTC - the same dates shown on DEA Maps are shown in local time (Sentinel-2 usually overpasses at around 11am in the morning).

Just so I'm understanding you correctly, you want to make your observations regularly spaced through time? So you get one timestep every 5 days even if the satellite data contains more frequent or less frequent images than that?

UQ-Raja-Ram-Aryal commented 7 months ago

I am working on monitoring and evaluating the Land Surface Phenology of Australian vegetation using Sentinel-2. So, the Time is more valuable for this work. In the same location how it will visit with two days. Those dates are in AM, why? So, How Can I adjust my data with a 5-day interval?

The times you get on data loaded from the datacube are in UTC - the same dates shown on DEA Maps are shown in local time (Sentinel-2 usually overpasses at around 11am in the morning).

Just so I'm understanding you correctly, you want to make your observations regularly spaced through time? So you get one timestep every 5 days even if the satellite data contains more frequent or less frequent images than that? image I am seeing Am here which is not supposed to be.

Is it good If I remove the intermediate one to maintain five-days interval data?

BexDunn commented 6 months ago

Hi @UQ-Raja-Ram-Aryal as @robbibt mentioned above, the times you are loading from the datacube are in UTC. The time displayed on DEA Maps is in local time, which I am guessing if you are at UQ will be UTC+10. Therefore Sentinel 2 overpasses will show up as a UTC time of roughly 1am, which corresponds to roughly 11am AEST (UTC+10).

It looks to me like you are getting data from Sentinel 2 A and Sentinel 2 B in your results:

# Load available data from Sentinel
ds = load_ard(dc=dc,
              products=['ga_s2am_ard_3', 'ga_s2bm_ard_3'],
              min_gooddata=0.90,
              **query)
print(ds) 

If you change your load to only get one of the products e.g.:

# Load available data from Sentinel 2A only
ds = load_ard(dc=dc,
              products=['ga_s2am_ard_3'],
              min_gooddata=0.90,
              **query)
print(ds) 

Then you may find that you have fewer results for the site as desired. Alternatively, you may be in the side lap of the swath and get results for the day and also for the following day.

BexDunn commented 5 months ago

I'm going to close this one as it isn't a bug. Happy to reopen if people think we need to improve our documentation or wish to discuss.