-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgrace_to_geotiff.py
154 lines (105 loc) · 4.61 KB
/
grace_to_geotiff.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
# -*- coding: utf-8 -*-
"""GRACE_to_GeoTiff.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/1dzAzUKYNvA48RiHNUjPuD9vqO8vy6_7b
"""
!pip install xarray
!pip install rioxarray
import xarray as xr
import rioxarray as rio
import rasterio as rrio
import numpy as np
from datetime import datetime, timedelta
"""## Testing a small NetCDF file for correct operation"""
!wget https://www.unidata.ucar.edu/software/netcdf/examples/sresa1b_ncar_ccsm3-example.nc
#Open the NetCDF
#Download the sample from https://www.unidata.ucar.edu/software/netcdf/examples/sresa1b_ncar_ccsm3-example.nc
ncfile = xr.open_dataset('sresa1b_ncar_ccsm3-example.nc')
#Extract the variable
pr = ncfile['pr']
#(Optional) convert longitude from (0-360) to (-180 to 180) (if required)
pr.coords['lon'] = (pr.coords['lon'] + 180) % 360 - 180
pr = pr.sortby(pr.lon)
#Define lat/long
pr = pr.rio.set_spatial_dims('lon', 'lat')
#Check for the CRS
pr.rio.crs
#(Optional) If your CRS is not discovered, you should be able to add it like so:
pr.rio.set_crs("epsg:4326")
#Saving the file
pr.rio.to_raster(r"sresa1b_ncar_ccsm3-exampleGeoTIFF.tif")
"""## Downloading the GRACE Time Series Dataset (also in NetCDF format)"""
!wget http://download.csr.utexas.edu/outgoing/grace/RL06_mascons/CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc
# this is a 783 MB file
#Open the NetCDF
rds = xr.open_dataset('CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc')
"""## Obtaining the data details and required parameters before outputting GeoTiff"""
interm_file_path = r"CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.tif"
final_file_path = r"CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02_2.tif"
print(rds)
#Extract the variable
pr = rds['lwe_thickness']
#(Optional) convert longitude from (0-360) to (-180 to 180) (if required)
pr.coords['lon'] = (pr.coords['lon'] + 180) % 360 - 180
pr = pr.sortby(pr.lon)
#Define lat/long
pr = pr.rio.set_spatial_dims('lon', 'lat')
#Check for the CRS
pr.rio.crs
#(Optional) If your CRS is not discovered, you should be able to add it like so:
pr.rio.write_crs("epsg:4326")
#Get the geotransform
transform = pr.rio.transform()
#Saving the file
pr.rio.to_raster(interm_file_path)
"""## Formatting the time information according to the start date"""
# Start date (in 'YYYY-MM-DD' format - verify correctness of this value for your dataset)
start_date = '2002-04-01'
# Gets the times to assign to each band later
times = pr.coords['time']
# Variable that stores the properly formatted datetimes
dts = []
start_offset = 0.0
start_dt = datetime.fromisoformat(start_date)
for i in range(len(times)):
decimal_offset = times[i].data.item()
current_dt = start_dt + timedelta(days=decimal_offset)
dts.append(datetime.isoformat(current_dt))
print(dts) # to verify that we have the right datetime format
"""## Renaming the bands according to time information; adding CRS"""
# Adding coordinate system information
crs = rrio.crs.CRS({"init": "epsg:4326"})
# Opening file to prepare to make a copy that we will add transform, crs information, and time information from earlier
src = rrio.open(interm_file_path, mode='r+')
# Opening the following will produce a warning about a lack of georeferencing, which will be resolved down below in write_transform.
with rrio.open(final_file_path, mode='w', driver='GTiff', width=src.width, height=src.height, count=src.count, dtype='float32') as dst:
for i in range(len(dts)): # goes through each band which is one per each time
# write the contents, plus the time as the description (name) of each band
band_data = src.read(i+1)
dst.write_band(i+1, band_data)
dst.set_band_description(i+1, dts[i])
# sets the CRS and transform
dst.crs = crs
dst.transform = transform
# closes the opened file
src.close()
"""## Verifying that the time series data are within the GeoTiff"""
from osgeo import gdal
import matplotlib.pyplot as plt
ds = gdal.Open(final_file_path)
ds_arr = ds.ReadAsArray()
print(ds_arr.shape)
print(ds.GetGeoTransform())
print(ds.GetProjection())
geoTransform = ds.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * ds.RasterXSize
miny = maxy + geoTransform[5] * ds.RasterYSize
print([minx, miny, maxx, maxy])
im = plt.imshow(ds_arr[0])
# this picture and the next one are flipped vertically because the geotransform has not been applied in plotting the raw array. It will appear properly as a GeoTiff.
im = plt.imshow(ds_arr[100])
# Enable faster download of file from Google Colab by zipping it first
!zip 'zippedtif.zip' 'CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02_2.tif'