2

For a given latitude/longitude location (for example, for a home or commercial property), I need to find out the storm surge at that location given a hurrican category.

The data is available in interactive form here: https://experience.arcgis.com/experience/203f772571cb48b1b8b50fdcc3272e2c/page/Category-4/

The same data is available as zipped GeoTIFF files at https://www.nhc.noaa.gov/nationalsurge/#data. There is one .tif file for each category, and several accompanying files for each.

I am able to load the data using the code below, but unable to do anything after as the array size is too large.
enter image description here Both data.plot() and data.values give me memory errors, even for the smallest of the files I get

"MemoryError: Unable to allocate 32.1 GiB for an array with shape (1, 136550, 252315) and data type uint8".

My laptop has 32GB RAM.

Is there a way to convert these files to GeoPandas format so I could simply use 'contains' to determine the property's expected storm surge?

5
  • 2
    You don't need to read the entire raster in, nor do you need to create a geodataframe from it. Just use the dataset.sample method. Here's an example geopandas.org/en/stable/gallery/geopandas_rasterio_sample.html
    – user2856
    Oct 8, 2022 at 21:15
  • Thank you for this - it works. Let me post final results once I am done. THanks a lot for helping a newcomer to GIS data.
    – Topchi
    Oct 9, 2022 at 23:22
  • When I use dataset.sample method, I need to use GPS coordinates per my TIFF file's CRS. However I am a returned a "no value" (255) unless the coordinates are exactly as coded in the file. For example: [int(x) for x in data.sample([(-74.018, 40.765)])] returns a value of 5, but the same for (-74.017, 40.765) returns a missing value equivalent 255. But these points are 100 feet away, and colored the same on the map equivalent. <br><br> Since I will provide any known coordinates, I would like the closest available point's value to be returned. Is there a way to do that?
    – Topchi
    Oct 10, 2022 at 0:03
  • 1
    Yes, the doc state "Values are from the nearest pixel. They are not interpolated." so if the point falls in a pixel, you should get a value. If not, can you ask a new question with code to demonstrate, example point coordinates and raster details (extent, rows, cols and CRS) so we can try to replicate.
    – user2856
    Oct 10, 2022 at 1:35
  • Thank you for this pointer as well. Yes, I am getting a value - there was a mistake in what I reported earlier (40.765, -74.017 was incorrect).
    – Topchi
    Oct 10, 2022 at 2:47

1 Answer 1

1

I was able to get to what I was trying to do. Many thanks to user2856 for the pointers. Here is the code I used

# Get lat and long from python using GeoPy
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="my-test-user-agent")
location = geolocator.geocode("33 Peck Slip, new york ny")

print(location.address)
print((location.longitude, location.latitude))

coords = [(location.longitude, location.latitude)]

# Use Cat1 SLOSH file from https://www.nhc.noaa.gov/nationalsurge/#data.
import rasterio
data = rasterio.open(r"us_Category1_MOM_Inundation_HIGH.tif")

print("\n Storm surge in feet:")
print([int(x) for x in data.sample(coords)])

enter image description here

1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.