Skip to content

Instantly share code, notes, and snippets.

@jmettes
Created February 4, 2020 00:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jmettes/96b64105fa9d0aced7321cef7d85748d to your computer and use it in GitHub Desktop.
Save jmettes/96b64105fa9d0aced7321cef7d85748d to your computer and use it in GitHub Desktop.
import cartopy.crs as ccrs
from mpl_toolkits.axes_grid1 import make_axes_locatable
proj = ccrs.PlateCarree()
fig, ax = plt.subplots(subplot_kw=dict(projection=proj), figsize=(10, 10))
s = ax.scatter(x['longitude'], x['latitude'], 10, x['iwv'])
gl = ax.gridlines(crs=proj, alpha=0.5, linestyle='--', draw_labels=True)
ax.coastlines(resolution='50m', color='black', linewidth=1)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05, axes_class=plt.Axes)
fig.colorbar(s, cax=cax)
import rasterio
from rasterio.plot import show
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
proj = ccrs.PlateCarree()
fig, ax = plt.subplots(subplot_kw=dict(projection=proj), figsize=(10, 10))
src = rasterio.open('20170208-20170220_VV_8rlks_eqa_small.unw.tif')
data = src.read()
data[data == 0] = np.nan
im = show(data, transform=src.transform, ax=ax, origin='upper')
cl = ax.coastlines()
cb = plt.colorbar(im.images[0])
plt.title('Interferogram)
plt.plot()
plt.savefig('interferogram.png')
gdal_translate -co COMPRESS=LZW -outsize 20% 20% 20170731-20170812_VV_8rlks_eqa.unw.tif 20170731-20170812_VV_8rlks_eqa_small.unw.tif
import xarray as xr
insar = xr.open_rasterio('../../../camden_data/geotiffs/S1/CAMDEN/T147D/20170731-20170812_VV_8rlks_eqa.unw.small.tif')
insar_point = insar.sel(y=latitude, x=longitude, method='nearest').item()
def bounding_box(lat, lon, offset_m=500):
earth = 6378.137 # radius of the earth in kilometer
m = (1 / ((2 * math.pi / 360) * earth)) / 1000 # 1 meter in degree
up = lat - (offset_m * m)
down = lat + (offset_m * m)
left = lon - (offset_m * m) / math.cos(lat * (math.pi / 180))
right = lon + (offset_m * m) / math.cos(lat * (math.pi / 180))
return up, down, left, right
up, down, left, right = bounding_box(lat=lat, lon=lon, offset_m500)
d = xr.open_rasterio('20151127-20151221_VV_8rlks_eqa.unw.tif')
region = d.sel(x=(d.x > left) & (d.x < right), y=(d.y > up) & (d.y < down))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment