topography extents

[Jonathan Burton]

Hi everyone,

I'm trying to use geoclaw to model an outburst flood. I have a high-res topography file for the drainage of an outburst flood and the DEM is just the immediate drainage valley (see screenshot attached). I input the full extent of the file in the setrun.py and when I run 

$ make .output

it says it reads the ascii .txt file but returns

 **** topo arrays do not cover domain

 **** area of overlap =    0.0000000000000000     

 **** area of domain  =   -4.8373999999998779E-002

Do I need a secondary DEM file to completely fill the designated extent or any ideas where this error is originating?

[Randall J LeVeque]

Hi Jonathan,

It is suspicious that it reports the area of your computational domain as negative: that should just be the area of the rectangle specified by clawdata.lower and clawdata.upper in setrun.py.   So that's the first thing to check, and make sure that the domain you specify is completely covered by the union of the topo files.

 - Randy

On Thu, Oct 19, 2023 at 10:55 AM Jonathan Burton <jonnybur...@gmail.com> wrote:

Hi everyone, I'm trying to use geoclaw to model an outburst flood. I have a high-res topography file for the drainage of an outburst flood and the DEM is just the immediate drainage valley (see screenshot attached). I input the full extent of

ZjQcmQRYFpfptBannerStart

This Message Is From an Untrusted Sender

You have not previously corresponded with this sender.

See https://itconnect.uw.edu/email-tags for additional information. Please contact the UW-IT Service Center, he...@uw.edu 206.221.5000, for assistance.

 

ZjQcmQRYFpfptBannerEnd

--
You received this message because you are subscribed to the Google Groups "claw-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to claw-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/claw-users/f4436b75-f6f3-4776-8e16-a54e146c68dcn%40googlegroups.com.

[Jonathan Burton]

The region is in southern and western hemisphere so the computational domain I set as:

#long

clawdata.lower[0] = -73.32

clawdata.upper[0] = -72.959

#Lat

clawdata.lower[1] = -46.899

clawdata.upper[1] = -47.033

These are the respective extents of the DEM file attached above, but everything in teal is noData, which may be my issue. I'll add another topofile to fill in the rest of the computational domain.

[Jon Elvar Wallevik]

OK, here is a silly question/comment:

In your setrun.py you are using latitude, longitude, i.e. coordinate_system = 2.
Is your map in perhaps another projection and you need to use coordinate_system = 1.
Check this on your original geodata (i.e. on the tif file, and not on the .asc file), with "gdalinfo filename.tif"

hope this helps
J

[Randall J LeVeque]

It looks like your lower latitude is north of your upper latitude, maybe these should be switched?  

Randy

Sent from Gmail Mobile

On Thu, Oct 19, 2023 at 7:12 PM Jonathan Burton <jonnybur...@gmail.com> wrote:

The region is in southern and western hemisphere so the computational domain I set as: #long clawdata. lower[0] = -73. 32 clawdata. upper[0] = -72. 959 #Lat clawdata. lower[1] = -46. 899 clawdata. upper[1] = -47. 033 These are the respective extents

To view this discussion on the web visit https://groups.google.com/d/msgid/claw-users/842cc9ad-1dff-467b-89dd-e045f12c416bn%40googlegroups.com.

[Jonathan Burton]

J - double checked on this and it is projected in WGS 1984, I switched the latitudes like Randy suggested (thank you, Randy, for pointing that out) and it is now reading in the topo file extents, but is all black. I think this may be an issue with the tif extent since I have an AOI within the extent boundaries and NoData to fill out the rest of the rectangle (see above). I'm going to try and overlay my higher-res tif with a low-res DEM to fill out the extents with data and see if that helps. Thanks all for the suggestions and help.

[Kyle Mandli]

I might suggest taking a look at the Topography class as it should be able to read in a topography file independent of what GeoClaw itself does.  In particular, it should be able to plot it.

import clawback.geoclaw.topo_tools as topo_tools
topo = topo_tools.Topography(path=PATH_TO_TOPO)
topo.plot()

This should also give you access to the information in the file.

The other question about the black seems odd to me, like you are trying to draw the grid cells in the plot and there are too many of them.  Also the background of the entire plot has a colorbar color so I am wondering if there may be something wrong with the plotting.

Kyle

On Oct 22, 2023 at 3:48 PM -0400, Jonathan Burton <jonnybur...@gmail.com>, wrote:

J - double checked on this and it is projected in WGS 1984, I switched the latitudes like Randy suggested (thank you, Randy, for pointing that out) and it is now reading in the topo file extents, but is all black. I think this may be an issue with the tif extent since I have an AOI within the extent boundaries and NoData to fill out the rest of the rectangle (see above). I'm going to try and overlay my higher-res tif with a low-res DEM to fill out the extents with data and see if that helps. Thanks all for the suggestions and help.

To view this discussion on the web visit https://groups.google.com/d/msgid/claw-users/ecfc7432-5513-4abe-84c6-994028b8b8f6n%40googlegroups.com.

[Jonathan Burton]

Hi Kyle, 

Thanks for the recommendations. I did pull the setrun.py from the Turzewski et al 2019 paper and amended it to my data, so that is a possibility with the grid cells since they look at a much larger region. I tried downscaling the number of grid cells in the setrun.py, but it's still not showing topo data, still just black. The command line prompt says it's reading in the two ascii topo files I have after 'make .output' command and both are ESRI ascii (type 3), but (full disclosure) I'm not the strongest at command-line code, trying to figure things out as I go. I created the setplot.py based off the Chile example and tried changing it for an outburst flood model, I just hadn't changed the colors yet (I have now).

[Kyle Mandli]

Hi Jonathan,

Here is some code from a very simple script if that would be helpful:

#!/usr/bin/env python

import os
import matplotlib.pyplot as plt
import clawpack.geoclaw.topotools as topotools

path = os.path.join(base_path, "ny_area.tt3")
topo_file = topotools.Topography(path, topo_type=3)
topo_file.read()

fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)

axes.pcolor(topo_file.X, topo_file.Y, topo_file.Z)

plt.show()

Of course change the path and topo_type to match what you have.  If you know a bit of matplotlib you can manipulate what it is that is being plotted.  This can be pasted into a script and run from the command line, from a Jupyter notebook, or from something like VSCode.  The snippet that I sent before is simpler of course.

If you want some eyes on the setplot.py I would say attach it to your next message and we can try and see if there's something glaring.

Kyle

To view this discussion on the web visit https://groups.google.com/d/msgid/claw-users/a00e3465-a1ae-42da-9d06-b18dba36dcc6n%40googlegroups.com.

[Jonathan Burton]

Hi Kyle,

I ran the code you gave me and it looks like it was able to read the topo ascii files so I don't think it's an inability to read and plot the topo file but it might be something related to what you originally suggested with the AMR grids because the plots are coming in now a bit funky, I attached a plot output for reference. There aren't any glaring warnings on the command code line when I'm running make commands, so it feels like I've just parameterized something incorrectly. Thanks for offering to take a loot at the setplot, that's where I feel least familiar.

[Kyle Mandli]

The biggest thing I think I saw was that plot_var = 0 corresponds to depth rather than sea-surface unless you have changed something.  Sea-surface is included in the output in variable 3 (zero-indexed) and is probably what you meant to have.  Also note that this variable is set to the topography if the cell is dry.  You can take a look at ways to handle this in $CLAW/geoclaw/src/python/geoclaw/surge/plot.py.

Kyle

To view this discussion on the web visit https://groups.google.com/d/msgid/claw-users/e253c54d-6f68-45f7-8e74-e922582d448fn%40googlegroups.com.