"""grdcut - Extract subregion from a grid or image or a slice from a cube."""fromtypingimportLiteralimportxarrayasxrfrompygmt.clibimportSessionfrompygmt.exceptionsimportGMTInvalidInputfrompygmt.helpersimport(build_arg_list,data_kind,fmt_docstring,kwargs_to_strings,use_alias,)__doctest_skip__=["grdcut"]
[docs]@fmt_docstring@use_alias(R="region",J="projection",N="extend",S="circ_subregion",V="verbose",Z="z_subregion",f="coltypes",)@kwargs_to_strings(R="sequence")defgrdcut(grid,kind:Literal["grid","image"]="grid",outgrid:str|None=None,**kwargs)->xr.DataArray|None:r""" Extract subregion from a grid or image or a slice from a cube. Produce a new ``outgrid`` file which is a subregion of ``grid``. The subregion is specified with ``region``; the specified range must not exceed the range of ``grid`` (but see ``extend``). If in doubt, run :func:`pygmt.grdinfo` to check range. Alternatively, define the subregion indirectly via a range check on the node values or via distances from a given point. Finally, you can give ``projection`` for oblique projections to determine the corresponding rectangular ``region`` that will give a grid that fully covers the oblique domain. Full option list at :gmt-docs:`grdcut.html` {aliases} Parameters ---------- {grid} kind The raster data kind. Valid values are ``"grid"`` and ``"image"``. When the input ``grid`` is a file name, it's difficult to determine if the file is a grid or an image, so we need to specify the raster kind explicitly. The default is ``"grid"``. {outgrid} {projection} {region} extend : bool or float Allow grid to be extended if new ``region`` exceeds existing boundaries. Give a value to initialize nodes outside current region. circ_subregion : str *lon/lat/radius*\[\ *unit*\][**+n**]. Specify an origin (*lon* and *lat*) and *radius*; append a distance *unit* and we determine the corresponding rectangular region so that all grid nodes on or inside the circle are contained in the subset. If **+n** is appended we set all nodes outside the circle to NaN. z_subregion : str [*min/max*\][**+n**\|\ **N**\|\ **r**]. Determine a new rectangular region so that all nodes outside this region are also outside the given z-range [-inf/+inf]. To indicate no limit on *min* or *max* only, specify a hyphen (-). Normally, any NaNs encountered are simply skipped and not considered in the range-decision. Append **+n** to consider a NaN to be outside the given z-range. This means the new subset will be NaN-free. Alternatively, append **+r** to consider NaNs to be within the data range. In this case we stop shrinking the boundaries once a NaN is found [Default simply skips NaNs when making the range decision]. Finally, if your core subset grid is surrounded by rows and/or columns that are all NaNs, append **+N** to strip off such columns before (optionally) considering the range of the core subset for further reduction of the area. {verbose} {coltypes} Returns ------- ret Return type depends on whether the ``outgrid`` parameter is set: - :class:`xarray.DataArray` if ``outgrid`` is not set - ``None`` if ``outgrid`` is set (grid output will be stored in the file set by ``outgrid``) Example ------- >>> import pygmt >>> # Load a grid of @earth_relief_30m data, with a longitude range of >>> # 10° E to 30° E, and a latitude range of 15° N to 25° N >>> grid = pygmt.datasets.load_earth_relief( ... resolution="30m", region=[10, 30, 15, 25] ... ) >>> # Create a new grid from an input grid, with a longitude range of >>> # 12° E to 15° E and a latitude range of 21° N to 24° N >>> new_grid = pygmt.grdcut(grid=grid, region=[12, 15, 21, 24]) """ifkindnotin{"grid","image"}:msg=f"Invalid raster kind: '{kind}'. Valid values are 'grid' and 'image'."raiseGMTInvalidInput(msg)# Determine the output data kind based on the input data kind.matchinkind:=data_kind(grid):case"grid"|"image":outkind=inkindcase"file":outkind=kindcase_:msg=f"Unsupported data type {type(grid)}."raiseGMTInvalidInput(msg)withSession()aslib:with(lib.virtualfile_in(check_kind="raster",data=grid)asvingrd,lib.virtualfile_out(kind=outkind,fname=outgrid)asvoutgrd,):kwargs["G"]=voutgrdlib.call_module(module="grdcut",args=build_arg_list(kwargs,infile=vingrd))returnlib.virtualfile_to_raster(vfname=voutgrd,kind=outkind,outgrid=outgrid)