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.

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

WKT and Shape files #22

Open
edgraffigna opened this issue Nov 27, 2024 · 5 comments
Open

WKT and Shape files #22

edgraffigna opened this issue Nov 27, 2024 · 5 comments

Comments

@edgraffigna
Copy link

Within the 3_Acces_HyP3_Data module, I've attempted to select the region of interest by creating a shapefile with polygon geometry using QGIS. I've generated shapefiles in both WGS84 and UTM Zone 19 South, but neither has yielded the desired results. I'm encountering an error. Could you please clarify the correct procedure for generating these files?
Specifically for the WKT file, I've copied the defined polygon directly from the ASF Vertex page and saved it as a .wkt text file. Furthermore, I've tried inputting this WKT into the request line, but both attempts have resulted in an error indicating that the specified area surpasses the dataset's boundaries.

@nshobert
Copy link

nshobert commented Feb 7, 2025

Hi @edgraffigna, did you ever get past this issue? I'm having trouble with the using the map and with the WKT approach. Can't get this subset step to work!

@Alex-Lewandowski
Copy link
Contributor

Hello,
@edgraffigna, I am sorry that I missed the issue you opened and left it hanging! @nshobert, thanks for giving this a bump!

The notebook expects a string WKT polygon or a shapefile. In either case, you will need to provide the coordinates in the projection of the data. All 4 shapefiles are needed. The notebook requests that the user selects the *.shp file in file browser; however, the accompanying *.dbf, *.prj, and *.shx files must also be present in the same directory.

There may be some occasions where you end up with a data stack in two projections. The second code cell in section 3 of the notebook checks and projects everything the predominantly represented projection. I would double-check that any shapefile or WKT you generate matches the data projection after this step.

Vertex AOI selection uses WGS84, so copy/pasting WKT directly to the notebook will not work. However, I understand that this is a natural flow for people working with ASF tools. I think it would be worth adding a step to the notebook to allow users to provide WKT in WGS84 or UTM.

@nshobert
Copy link

nshobert commented Feb 7, 2025

Thanks for the quick response, @Alex-Lewandowski. I'm trying the shapefile workflow and get NameError: name 'epsg' is not defined. Is it possible that espg is supposed to be outside of this else statement? 2nd to last line

else:
    correct_wkt_input = False
    while not correct_wkt_input:
        if common_coverage:
            wkt = (f'POLYGON(({common_extents[0]} {common_extents[1]}, {common_extents[2]} {common_extents[1]}, {common_extents[2]} '
                   f'{common_extents[3]}, {common_extents[0]} {common_extents[3]}, {common_extents[0]} {common_extents[1]}))')
            wkt_shapely_geom = shapely.wkt.loads(wkt)
        else:
            wkt, wkt_shapely_geom = util.get_valid_wkt()
            
        wkt_ogr_geom = ogr.CreateGeometryFromWkt(wkt)
        
        if not util.check_within_bounds(wkt_shapely_geom, gdf):
            print('WKT exceeds bounds of at least one dataset')
            if common_coverage:
                raise Exception('Error determining area of common coverage')
            else:
                continue
    
        correct_wkt_input = True
        
    shp_path = data_path / f'shape_{datetime.strftime(datetime.now(), "%Y%m%dT%H%M%S")}.shp'
    epsg = int(gdf.iloc[0]['EPSG'])
    util.save_shapefile(wkt_ogr_geom, epsg, shp_path)

@Alex-Lewandowski
Copy link
Contributor

Hi @nshobert, In this code cell, the shapefile workflow would take the elif branch above the else branch that you showed. In that elif, you will be prompted to select a shapefile from a file browser:

elif shapefile:
    print('Select a shapefile (*.shp)')
    shp_fc = FileChooser(Path.home())
    display(shp_fc)

The else branch is for the WKT workflow, which includes a step that produces a shapefile used for subsetting. The epsg variable should remain in that branch. It should only create a new shapefile if you are using the WKT option.

The else (WKT) branch:

  1. creates a shapely geometry from the WKT
  2. validates that the shape falls within the bounds of the data
  3. generates the shapefile from the geometry, which is used for subsetting 2 code cells later.

@nshobert
Copy link

nshobert commented Feb 9, 2025

Thanks, @Alex-Lewandowski. How is epsg assigned if we don't enter the else (WKT) branch? The epsg is called in the next cell in if not draw branch, which is true for the user-provided shapefile workflow:

plt.close()
if shapefile:
    shp_path = Path(shp_fc.selected)
    if shp_path.suffix != '.shp':
        raise Exception(f'Selected file suffix not ".shp"')
if not draw:   
    shp_gdf = gpd.read_file(shp_path)
    shp_gdf = shp_gdf.to_crs(crs='EPSG:4326')

    box1 = gpd.GeoDataFrame({"geometry": [box(*max_extents)]}, crs=f'EPSG:{epsg}')

I get NameError: name 'epsg' is not defined

I hope I'm not missing something obvious 😅

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants