Skip to content

Commit

Permalink
debug new epsg questionnaire
Browse files Browse the repository at this point in the history
  • Loading branch information
@BuddyVolly
BuddyVolly committed Sep 8, 2021
1 parent 6383d52 commit f9e629e
Show file tree
Hide file tree
Showing 5 changed files with 85 additions and 153 deletions.
54 changes: 34 additions & 20 deletions ost/Project.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -662,7 +662,9 @@ def bursts_to_ards(
self.ard_parameters['single_ARD']['dem'][
'dem_name'] = 'Copernicus 30m Global DEM'

if self.ard_parameters['dem']['out_projection'] == 4326:
if self.ard_parameters['single_ARD']['dem'][
'out_projection'
] == 4326:

logger.info(
'The scene\'s location is towards the poles. '
Expand Down Expand Up @@ -819,26 +821,38 @@ def grds_to_ards(
# when outside SRTM
center_lat = loads(self.aoi).centroid.y
if float(center_lat) > 59 or float(center_lat) < -59:
logger.info(
'Scene is outside SRTM coverage. Snap will therefore use '
'the GETASSE30 DEM. Also consider to use a stereographic '
'projection. and project to NSIDC Sea Ice Polar '
'Stereographic North projection (EPSG 3413).'
)
epsg = input(
'Please type the EPSG you want to project the output data or '
'just press enter for keeping Lat/Lon coordinate system '
'(e.g. 3413 for NSIDC Sea Ice Polar Stereographic North '
'projection, or 3976 for NSIDC Sea Ice Polar Stereographic '
'South projection'
)
if not epsg:
epsg = 4326

self.ard_parameters['single_ARD']['dem'][
'dem_name'] = 'GETASSE30'
self.ard_parameters['single_ARD']['dem'][
'out_projection'] = int(epsg)
if 'SRTM' in self.ard_parameters['single_ARD']['dem']['dem_name']:
logger.info(
'Scene is outside SRTM coverage. Snap will therefore use '
'the Copernicus 30m Global DEM. '
)

self.ard_parameters['single_ARD']['dem'][
'dem_name'] = 'Copernicus 30m Global DEM'

if self.ard_parameters['single_ARD']['dem'][
'out_projection'
] == 4326:

logger.info(
'The scene\'s location is towards the poles. '
'Consider to use a stereographic projection.'
)

epsg = input(
'Type an alternative EPSG code for the projection of the '
'output data or just press enter for keeping Lat/Lon '
'coordinate system (e.g. 3413 for NSIDC Sea Ice Polar '
'Stereographic North projection, or 3976 for '
'NSIDC Sea Ice Polar Stereographic South projection'
)

if not epsg:
epsg = 4326

self.ard_parameters['single_ARD']['dem'][
'out_projection'] = int(epsg)

# --------------------------------------------
# 3 subset determination
Expand Down
153 changes: 30 additions & 123 deletions ost/generic/ts_ls_mask.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,122 +5,23 @@
import shutil
import logging
from pathlib import Path
from tempfile import TemporaryDirectory

import rasterio
import numpy as np
import geopandas as gpd
from shapely.ops import unary_union
from retrying import retry

try:
import gdal
except ModuleNotFoundError as e:
from osgeo import gdal
except ModuleNotFoundError as e:
raise e
try:
from osgeo import gdal
except ModuleNotFoundError:
raise e

from ost.helpers import raster as ras, vector as vec

logger = logging.getLogger(__name__)

@retry(stop_max_attempt_number=3, wait_fixed=1)
def mt_layover_old(list_of_files, config_file):
"""
:param list_of_files:
:param config_file:
:return:
"""

# this is a godale thing
with open(config_file) as file:
config_dict = json.load(file)
temp_dir = Path(config_dict['temp_dir'])
update_extent = (
config_dict['processing']['time-series_ARD']['apply_ls_mask']
)

target_dir = Path(list_of_files[0]).parent.parent.parent
outfile = target_dir.joinpath(f'{target_dir.name}.ls_mask.tif')
extent = target_dir.joinpath(f'{target_dir.name}.extent.gpkg')
burst_dir = Path(outfile).parent
burst = burst_dir.name

logger.info(
f'Creating common Layover/Shadow mask for track {target_dir.name}.'
)

with TemporaryDirectory(prefix=f'{temp_dir}/') as temp:

# temp to Path object
temp = Path(temp)

# create path to temp file
ls_layer = temp.joinpath(Path(outfile).name)

# create a vrt-stack out of
vrt_options = gdal.BuildVRTOptions(srcNodata=0, separate=True)
gdal.BuildVRT(
str(temp.joinpath('ls.vrt')),
list_of_files,
options=vrt_options
)

with rasterio.open(temp.joinpath('ls.vrt')) as src:

# get metadata
meta = src.meta
# update driver and reduced band count
meta.update(driver='GTiff', count=1, dtype='uint8')

# create outfiles
with rasterio.open(ls_layer, 'w', **meta) as out_min:

# loop through blocks
for _, window in src.block_windows(1):

# read array with all bands
stack = src.read(range(1, src.count + 1), window=window)

# get stats
arr_max = np.nanmax(stack, axis=0)
arr = np.divide(arr_max, arr_max)

out_min.write(np.uint8(arr), window=window, indexes=1)

ras.mask_by_shape(
ls_layer, outfile, extent, to_db=False,
datatype='uint8', rescale=False, ndv=0
)

ls_layer.unlink()

extent_ls_masked = None
if update_extent:

logger.info(
'Calculating symmetrical difference of extent and ls_mask'
)

# polygonize the multi-temporal ls mask
ras.polygonize_raster(outfile, f'{str(outfile)[:-4]}.gpkg')

# create file for masked extent
extent_ls_masked = burst_dir.joinpath(
f'{burst}.extent.masked.gpkg'
)

# calculate difference between burst extent
# and ls mask, for masked extent
try:
vec.difference(
extent, f'{outfile.stem}.gpkg', extent_ls_masked
)
except:
shutil.copy(extent, extent_ls_masked)

return burst_dir, list_of_files, outfile, extent_ls_masked

@retry(stop_max_attempt_number=3, wait_fixed=1)
def mt_layover(list_of_ls):
Expand All @@ -137,34 +38,40 @@ def mt_layover(list_of_ls):
f'Creating common Layover/Shadow mask for track {target_dir.name}.'
)

y = 0
for i, file in enumerate(list_of_ls):

if i == 0:
if y == 0:
df1 = gpd.read_file(file)
df1 = df1[~(df1.geometry.is_empty | df1.geometry.isna())]
geom = df1.geometry.buffer(0).unary_union
if len(df1) > 0:
geom = df1.geometry.buffer(0).unary_union
y = 1

if i > 0:
if y > 0:

df2 = gpd.read_file(file)
df2 = df2[~(df2.geometry.is_empty | df2.geometry.isna())]
geom2 = df2.geometry.buffer(0).unary_union
geom = unary_union([geom, geom2])

# make geometry valid in case it isn't
geom = geom.buffer(0)

# remove slivers
buffer = 0.00001
geom = geom.buffer(
-buffer, 1, join_style=2
).buffer(
buffer, 1, cap_style=1, join_style=2
).__geo_interface__

# write to output
with open(outfile, 'w') as file:
json.dump(geom, file)

# create difference file for valid data shape
vec.difference(bounds, outfile, valid_file)
if y > 0:
# make geometry valid in case it isn't
geom = geom.buffer(0)

# remove slivers
buffer = 0.00001
geom = geom.buffer(
-buffer, 1, join_style=2
).buffer(
buffer, 1, cap_style=1, join_style=2
).__geo_interface__

# write to output
with open(outfile, 'w') as file:
json.dump(geom, file)

# create difference file for valid data shape
vec.difference(bounds, outfile, valid_file)
else:
shutil.copy(bounds, valid_file)
18 changes: 11 additions & 7 deletions ost/helpers/vector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,10 @@
try:
import ogr
except ModuleNotFoundError as e:
from osgeo import ogr
except ModuleNotFoundError as e:
raise e
try:
from osgeo import ogr
except ModuleNotFoundError:
raise e

from shapely.ops import transform
from shapely.wkt import loads
Expand Down Expand Up @@ -61,9 +62,9 @@ def aoi_to_wkt(aoi):
if Path(aoi).exists():
try:
gdf = gpd.GeoDataFrame.from_file(aoi)
if gdf.crs != {'epsg:4326'}:
if gdf.crs != 'epsg:4326':
try:
gdf = gdf.geometry.to_crs({'epsg:4326'})
gdf = gdf.geometry.to_crs('epsg:4326')
except:
raise ValueError('No valid OST AOI definition.')
# return AOI as single vector object
Expand Down Expand Up @@ -339,8 +340,6 @@ def wkt_to_gdf(wkt):
:return:
"""

warnings.filterwarnings('ignore', r'syntax is deprecated', FutureWarning)

# load wkt
geometry = loads(wkt)

Expand Down Expand Up @@ -427,6 +426,11 @@ def exterior(infile, outfile, buffer=None):

def difference(infile1, infile2, outfile):

import warnings
warnings.filterwarnings(
'ignore', 'Geometry is in a geographic CRS', UserWarning
)

gdf1 = gpd.read_file(infile1)
gdf2 = gpd.read_file(infile2)

Expand Down
7 changes: 6 additions & 1 deletion ost/s1/refine_inventory.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -490,9 +490,14 @@ def search_refinement(
"""

# creat AOI GeoDataframe and calulate area
warnings.filterwarnings(
'ignore', 'Geometry is in a geographic CRS', UserWarning
)

# create AOI GeoDataframe and calulate area
aoi_gdf = vec.wkt_to_gdf(aoi)
aoi_area = aoi_gdf.area.sum()

# get all polarisations apparent in the inventory
pols = inventory_df['polarisationmode'].unique()

Expand Down
6 changes: 4 additions & 2 deletions ost/s1/s1scene.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -717,7 +717,7 @@ def set_external_dem(self, dem_file, ellipsoid_correction=True):
import rasterio

# check if file exists
if not Path(dem_file).eixtst():
if not Path(dem_file).exists():
raise FileNotFoundError(f'No file found at {dem_file}.')

# get no data value
Expand Down Expand Up @@ -799,7 +799,9 @@ def create_ard(self, infile, out_dir, subset=None, overwrite=False):
self.ard_parameters['single_ARD']['dem'][
'dem_name'] = 'Copernicus 30m Global DEM'

if self.ard_parameters['dem']['out_projection'] == 4326:
if self.ard_parameters['single_ARD']['dem'][
'out_projection'
] == 4326:

logger.info(
'The scene\'s location is towards the poles. '
Expand Down

0 comments on commit f9e629e

Please sign in to comment.