Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: use clongdouble instead of longcomplex for numpy 2.0 #45

Merged
merged 2 commits into from
Jun 27, 2024
Merged

fix: use clongdouble instead of longcomplex for numpy 2.0 #45

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3f0f15d
fix: use `clongdouble` instead of `longcomplex` for numpy 2.0
tsutterley Jun 27, 2024
89dad11
feat: added option to not run with crossover measurements
tsutterley Jun 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions DEM/MPI_DEM_ICESat2_ATL03.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,8 @@
h5py = gz.utilities.import_dependency('h5py')
MPI = gz.utilities.import_dependency('mpi4py.MPI')
pyproj = gz.utilities.import_dependency('pyproj')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')
timescale = gz.utilities.import_dependency('timescale')

# digital elevation models
Expand Down Expand Up @@ -277,7 +278,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -493,7 +494,7 @@ def main():
X,Y = transformer.transform(longitude, latitude)

# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(list(zip(X[ind], Y[ind])))
xy_point = shapely.geometry.MultiPoint(list(zip(X[ind], Y[ind])))

# create complete masks for each DEM tile
associated_map = {}
Expand Down
7 changes: 4 additions & 3 deletions DEM/MPI_DEM_ICESat2_ATL06.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,8 @@
h5py = gz.utilities.import_dependency('h5py')
MPI = gz.utilities.import_dependency('mpi4py.MPI')
pyproj = gz.utilities.import_dependency('pyproj')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')
timescale = gz.utilities.import_dependency('timescale')

# digital elevation models
Expand Down Expand Up @@ -277,7 +278,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -493,7 +494,7 @@ def main():
X,Y = transformer.transform(longitude, latitude)

# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(np.c_[X[ind], Y[ind]])
xy_point = shapely.geometry.MultiPoint(np.c_[X[ind], Y[ind]])

# create complete masks for each DEM tile
associated_map = {}
Expand Down
7 changes: 4 additions & 3 deletions DEM/MPI_DEM_ICESat2_ATL11.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,8 @@
h5py = gz.utilities.import_dependency('h5py')
MPI = gz.utilities.import_dependency('mpi4py.MPI')
pyproj = gz.utilities.import_dependency('pyproj')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')
timescale = gz.utilities.import_dependency('timescale')

# digital elevation models
Expand Down Expand Up @@ -257,7 +258,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -470,7 +471,7 @@ def main():
fileID[ptx]['latitude'][:])

# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(np.c_[X[ind], Y[ind]])
xy_point = shapely.geometry.MultiPoint(np.c_[X[ind], Y[ind]])

# create complete masks for each DEM tile
associated_map = {}
Expand Down
7 changes: 4 additions & 3 deletions DEM/MPI_DEM_ICESat_GLA12.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,8 @@
MPI = gz.utilities.import_dependency('mpi4py.MPI')
pyproj = gz.utilities.import_dependency('pyproj')
pyTMD = gz.utilities.import_dependency('pyTMD')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')

# digital elevation models
elevation_dir = {}
Expand Down Expand Up @@ -238,7 +239,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -463,7 +464,7 @@ def main():
X,Y = transformer.transform(lon_40HZ, lat_40HZ)

# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(np.c_[X[ind], Y[ind]])
xy_point = shapely.geometry.MultiPoint(np.c_[X[ind], Y[ind]])

# create complete masks for each DEM tile
associated_map = {}
Expand Down
7 changes: 4 additions & 3 deletions DEM/MPI_interpolate_DEM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,8 @@
gdal = gz.utilities.import_dependency('osgeo.gdal')
pyproj = gz.utilities.import_dependency('pyproj')
pyTMD = gz.utilities.import_dependency('pyTMD')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')
timescale = gz.utilities.import_dependency('timescale')

# digital elevation models
Expand Down Expand Up @@ -292,7 +293,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -509,7 +510,7 @@ def main():
dem_h.mask = np.ones((n_pts),dtype=bool)

# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(np.c_[X[ind], Y[ind]])
xy_point = shapely.geometry.MultiPoint(np.c_[X[ind], Y[ind]])

# create complete masks for each DEM tile
associated_map = {}
Expand Down
7 changes: 4 additions & 3 deletions DEM/check_DEM_ICESat2_ATL06.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,8 @@
fiona = gz.utilities.import_dependency('fiona')
is2tk = gz.utilities.import_dependency('icesat2_toolkit')
pyproj = gz.utilities.import_dependency('pyproj')
geometry = gz.utilities.import_dependency('shapely.geometry')
shapely = gz.utilities.import_dependency('shapely')
shapely.geometry = gz.utilities.import_dependency('shapely.geometry')

# digital elevation models
elevation_dir = {}
Expand Down Expand Up @@ -185,7 +186,7 @@ def read_DEM_index(index_file, DEM_MODEL):
# extract Polar Stereographic coordinates for entity
x = [ul[0],ur[0],lr[0],ll[0],ul2[0]]
y = [ul[1],ur[1],lr[1],ll[1],ul2[1]]
poly_obj = geometry.Polygon(np.c_[x, y])
poly_obj = shapely.geometry.Polygon(np.c_[x, y])
# Valid Polygon may not possess overlapping exterior or interior rings
if (not poly_obj.is_valid):
poly_obj = poly_obj.buffer(0)
Expand Down Expand Up @@ -258,7 +259,7 @@ def check_DEM_ICESat2_ATL06(INPUT_FILE,
# convert projection from latitude/longitude to tile EPSG
X,Y = transformer.transform(longitude, latitude)
# convert reduced x and y to shapely multipoint object
xy_point = geometry.MultiPoint(np.c_[X, Y])
xy_point = shapely.geometry.MultiPoint(np.c_[X, Y])

# create complete masks for each DEM tile
intersection_map = {}
Expand Down
49 changes: 25 additions & 24 deletions DEM/create_GIMP_tile_index.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,9 +71,10 @@
import grounding_zones as gz

# attempt imports
gdal = gz.utilities.import_dependency('osgeo.gdal')
osr = gz.utilities.import_dependency('osgeo.osr')
ogr = gz.utilities.import_dependency('osgeo.ogr')
osgeo = gz.utilities.import_dependency('osgeo')
osgeo.gdal = gz.utilities.import_dependency('osgeo.gdal')
osgeo.osr = gz.utilities.import_dependency('osgeo.osr')
osgeo.ogr = gz.utilities.import_dependency('osgeo.ogr')

# PURPOSE: read GIMP image mosaic and create tile shapefile
def create_GIMP_tile_index(base_dir, VERSION, MODE=0o775):
Expand All @@ -94,30 +95,30 @@ def create_GIMP_tile_index(base_dir, VERSION, MODE=0o775):
parser = lxml.etree.HTMLParser()

# save DEM tile outlines to ESRI shapefile
driver = ogr.GetDriverByName('Esri Shapefile')
driver = osgeo.ogr.GetDriverByName('Esri Shapefile')
output_shapefile = ddir.joinpath(ff.format(VERSION,'shp'))
ds = driver.CreateDataSource(str(output_shapefile))
# set the spatial reference info
# EPSG: 3413 (NSIDC Sea Ice Polar Stereographic North, WGS84)
SpatialReference = osr.SpatialReference()
SpatialReference = osgeo.osr.SpatialReference()
SpatialReference.ImportFromEPSG(3413)
layer = ds.CreateLayer('', SpatialReference, ogr.wkbPolygon)
layer = ds.CreateLayer('', SpatialReference, osgeo.ogr.wkbPolygon)
# Add shapefile attributes (following attributes from ArcticDEM and REMA)
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
layer.CreateField(ogr.FieldDefn('name', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('tile', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('nd_value', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('resolution', ogr.OFTInteger))
layer.CreateField(ogr.FieldDefn('lastmodtm', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('fileurl', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('rel_ver', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('spec_type', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('reg_src', ogr.OFTString))
field_area = ogr.FieldDefn('st_area_sh', ogr.OFTReal)
layer.CreateField(osgeo.ogr.FieldDefn('id', osgeo.ogr.OFTInteger))
layer.CreateField(osgeo.ogr.FieldDefn('name', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('tile', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('nd_value', osgeo.ogr.OFTReal))
layer.CreateField(osgeo.ogr.FieldDefn('resolution', osgeo.ogr.OFTInteger))
layer.CreateField(osgeo.ogr.FieldDefn('lastmodtm', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('fileurl', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('rel_ver', osgeo.ogr.OFTReal))
layer.CreateField(osgeo.ogr.FieldDefn('spec_type', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('reg_src', osgeo.ogr.OFTString))
field_area = osgeo.ogr.FieldDefn('st_area_sh', osgeo.ogr.OFTReal)
field_area.SetWidth(24)
field_area.SetPrecision(10)
layer.CreateField(field_area)
layer.CreateField(ogr.FieldDefn('st_length_', ogr.OFTInteger))
layer.CreateField(osgeo.ogr.FieldDefn('st_length_', osgeo.ogr.OFTInteger))
defn = layer.GetLayerDefn()

# create a counter for shapefile id
Expand Down Expand Up @@ -147,8 +148,8 @@ def create_GIMP_tile_index(base_dir, VERSION, MODE=0o775):

# use GDAL memory-mapped file to read dem
mmap_name = f'/vsimem/{uuid.uuid4().hex}'
gdal.FileFromMemBuffer(mmap_name, fileID.read())
dataset = gdal.Open(mmap_name)
osgeo.gdal.FileFromMemBuffer(mmap_name, fileID.read())
dataset = osgeo.gdal.Open(mmap_name)

# get dimensions of tile
xsize = dataset.RasterXSize
Expand All @@ -172,10 +173,10 @@ def create_GIMP_tile_index(base_dir, VERSION, MODE=0o775):
st_length_ = int(2*(xmax-xmin) + 2*(ymax-ymin))
# close the dataset
dataset = None
gdal.Unlink(mmap_name)
osgeo.gdal.Unlink(mmap_name)

# Create a new feature (attribute and geometry)
feature = ogr.Feature(defn)
feature = osgeo.ogr.Feature(defn)
# Add shapefile attributes
feature.SetField('id', id)
feature.SetField('name', colname.replace('.tif',''))
Expand All @@ -192,11 +193,11 @@ def create_GIMP_tile_index(base_dir, VERSION, MODE=0o775):
feature.SetField('st_length_', st_length_)

# create LineString object and add x/y points
ring_obj = ogr.Geometry(ogr.wkbLinearRing)
ring_obj = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
for x, y in zip(xbox, ybox):
ring_obj.AddPoint(x, y)
# create Polygon object for LineString of tile
poly_obj = ogr.Geometry(ogr.wkbPolygon)
poly_obj = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
poly_obj.AddGeometry(ring_obj)
feature.SetGeometry(poly_obj)
layer.CreateFeature(feature)
Expand Down
53 changes: 27 additions & 26 deletions DEM/nsidc_convert_GIMP_DEM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,9 +72,10 @@
import grounding_zones as gz

# attempt imports
gdal = gz.utilities.import_dependency('osgeo.gdal')
osr = gz.utilities.import_dependency('osgeo.osr')
ogr = gz.utilities.import_dependency('osgeo.ogr')
osgeo = gz.utilities.import_dependency('osgeo')
osgeo.gdal = gz.utilities.import_dependency('osgeo.gdal')
osgeo.osr = gz.utilities.import_dependency('osgeo.osr')
osgeo.ogr = gz.utilities.import_dependency('osgeo.ogr')

# PURPOSE: read GIMP image mosaic and output as gzipped tar file
def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
Expand Down Expand Up @@ -146,8 +147,8 @@ def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
fileID.seek(0)
# use GDAL memory-mapped file to read dem
mmap_name = f"/vsimem/{uuid.uuid4().hex}"
gdal.FileFromMemBuffer(mmap_name, fileID.read())
dataset = gdal.Open(mmap_name)
osgeo.gdal.FileFromMemBuffer(mmap_name, fileID.read())
dataset = osgeo.gdal.Open(mmap_name)

# get dimensions of tile
xsize = dataset.RasterXSize
Expand All @@ -171,10 +172,10 @@ def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
st_length_ = int(2*(xmax-xmin) + 2*(ymax-ymin))
# close the dataset
dataset = None
gdal.Unlink(mmap_name)
osgeo.gdal.Unlink(mmap_name)

# save DEM tile outlines to ESRI shapefile
driver = ogr.GetDriverByName('Esri Shapefile')
driver = osgeo.ogr.GetDriverByName('Esri Shapefile')
# use GDAL memory-mapped file
mmap = {}
for key in ('dbf','prj','shp','shx'):
Expand All @@ -183,25 +184,25 @@ def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
ds = driver.CreateDataSource(mmap['shp'])
# set the spatial reference info
# EPSG: 3413 (NSIDC Sea Ice Polar Stereographic North, WGS84)
SpatialReference = osr.SpatialReference()
SpatialReference = osgeo.osr.SpatialReference()
SpatialReference.ImportFromEPSG(3413)
layer = ds.CreateLayer('', SpatialReference, ogr.wkbPolygon)
layer = ds.CreateLayer('', SpatialReference, osgeo.ogr.wkbPolygon)
# Add shapefile attributes (following attributes from ArcticDEM and REMA)
layer.CreateField(ogr.FieldDefn('DEM_ID', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('DEM_NAME', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('TILE', ogr.OFTString))
layer.CreateField(ogr.FieldDefn('ND_VALUE', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('DEM_RES', ogr.OFTInteger))
layer.CreateField(ogr.FieldDefn('REL_VER', ogr.OFTReal))
layer.CreateField(ogr.FieldDefn('REG_SRC', ogr.OFTString))
field_area = ogr.FieldDefn('ST_AREA', ogr.OFTReal)
layer.CreateField(osgeo.ogr.FieldDefn('DEM_ID', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('DEM_NAME', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('TILE', osgeo.ogr.OFTString))
layer.CreateField(osgeo.ogr.FieldDefn('ND_VALUE', osgeo.ogr.OFTReal))
layer.CreateField(osgeo.ogr.FieldDefn('DEM_RES', osgeo.ogr.OFTInteger))
layer.CreateField(osgeo.ogr.FieldDefn('REL_VER', osgeo.ogr.OFTReal))
layer.CreateField(osgeo.ogr.FieldDefn('REG_SRC', osgeo.ogr.OFTString))
field_area = osgeo.ogr.FieldDefn('ST_AREA', osgeo.ogr.OFTReal)
field_area.SetWidth(24)
field_area.SetPrecision(10)
layer.CreateField(field_area)
layer.CreateField(ogr.FieldDefn('ST_LENGTH', ogr.OFTInteger))
layer.CreateField(osgeo.ogr.FieldDefn('ST_LENGTH', osgeo.ogr.OFTInteger))
defn = layer.GetLayerDefn()
# Create a new feature (attribute and geometry)
feature = ogr.Feature(defn)
feature = osgeo.ogr.Feature(defn)
# Add shapefile attributes
feature.SetField('DEM_ID', f'{tile}_{res}')
feature.SetField('DEM_NAME', colname)
Expand All @@ -213,11 +214,11 @@ def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
feature.SetField('ST_AREA', st_area_sh)
feature.SetField('ST_LENGTH', st_length_)
# create LineString object and add x/y points
ring_obj = ogr.Geometry(ogr.wkbLinearRing)
ring_obj = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
for x,y in zip(xbox,ybox):
ring_obj.AddPoint(x,y)
# create Polygon object for LineString of tile
poly_obj = ogr.Geometry(ogr.wkbPolygon)
poly_obj = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
poly_obj.AddGeometry(ring_obj)
feature.SetGeometry(poly_obj)
layer.CreateFeature(feature)
Expand All @@ -234,16 +235,16 @@ def nsidc_convert_GIMP_DEM(base_dir, VERSION, MODE=0o775):
for key in ('dbf','prj','shp','shx'):
output = f'{tile}_{res}_index.{key}'
info4 = tarfile.TarInfo(name=posixpath.join(subdir,'index',output))
info4.size = gdal.VSIStatL(mmap[key]).size
info4.size = osgeo.gdal.VSIStatL(mmap[key]).size
info4.mtime = remote_mtime
info4.mode = MODE
fp = gdal.VSIFOpenL(mmap[key],'rb')
fp = osgeo.gdal.VSIFOpenL(mmap[key],'rb')
with io.BytesIO() as fileID:
fileID.write(gdal.VSIFReadL(1, info4.size, fp))
fileID.write(osgeo.gdal.VSIFReadL(1, info4.size, fp))
fileID.seek(0)
tar.addfile(tarinfo=info4,fileobj=fileID)
gdal.VSIFCloseL(fp)
gdal.Unlink(mmap[key])
osgeo.gdal.VSIFCloseL(fp)
osgeo.gdal.Unlink(mmap[key])

# close tar file
tar.close()
Expand Down
Loading