Hot answers tagged

5

You probably have a null geometry, try filtering them out. import geopandas as gpd import pandas as pd from rasterstats import zonal_stats shp = 'vector.shp' ras = 'raster.tif' gdf = gpd.read_file(shp) # filter out empty or null geometries gdf = gdf[~(gdf['geometry'].is_empty | gdf['geometry'].isna())] # zonal stats stats = zonal_stats(gdf['geometry'], ...


3

Your code works fine. If you look at SVDNB_npp_d20120301.vcld.tif (src.meta) you'll see nodata = 255 and you specified filled=True which means fill out_image with nodata values. So when you print your array, you are seeing the nodata values around the masked area of actual values. If you write out the image and load it in a GIS, you'll see: # To mask the ...


3

You can use pyqgis. I have a field called kkod and I'm adding a new field called new_kkod and calculating. Adjust layername, field names and type of new field, and the d dictionary of old and new values. You can add all the values you like to it. layer = QgsProject.instance().mapLayersByName('buildings')[0] #Add field pr = layer.dataProvider() pr....


3

Use Capa.getFeatures() in each loop without assigning a variable(Ejes_duplicados) contador = 0 for eje in Capa.getFeatures(): # <<< contador = contador+1 print("----------contador: ", contador) LoteID = eje['ID'] print("1_ Lote id: ", LoteID) Clase_Lote = eje['LOTE_CLASS'] print("2_ Clase Lote: ",...


3

If you do not have ArcGIS Pro currently installed on your machine If you installed the arcgis library through conda but do not have ArcGIS Pro installed on yout machine, you need to call the activate.bat file from the main conda folder and point it to your newly created environment. Edit your .bat file to look like this: @echo on set env_name=ArcGIS_Test ...


3

Joseph's answer is correct, although there needs to be another parameter included: "PREDICATE". ie: processing.run("qgis:selectbylocation", { "INPUT":lyr_input,\ "PREDICATE":0,\ "INTERSECT":lyr_intersect,\ "METHOD":0,\ "OUTPUT":path} ) From processing.algorithmHelp(&...


2

Here's a very naive way using rasterio.shutil.copy to generate VRT XML for each source raster, combine them and output a stacked VRT (equivalent to gdalbuildvrt -separate). By naive, I mean it does not check any of the properties and assumes rasters are all the same shape, extent, pixel size, data type, projection etc. import xml.etree.ElementTree as ET ...


2

Maybe I misunderstand you, but you can use .isin like this to create a third dataframe of the records in df1 that exists in df2: import pandas as pd data1 = {'col_1': [3, 2, 1, 0], 'col_2': ['a', 'b', 'c', 'd']} data2 = {'col_1': [9, 2, 1, 8], 'col_2': ['e', 'f', 'g', 'h']} df1 = pd.DataFrame.from_dict(data1) # col_1 col_2 #0 3 a #1 2 b ...


2

The answer from user30184 was correct. Adding target.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) after importing the SR from an EPSG code gave me the desired output! Here is the example code I posted before, showing how it is used. from osgeo import ogr,osr source = osr.SpatialReference() source.ImportFromEPSG(26913) target = osr.SpatialReference()...


1

Your polygon intersects with several different Sentinel 1 tiles. The UTM zone (shown in red) is probably taken from the center of a tile, in which case it makes sense they are different for different dates. You could reproject your images to WGS before exporting but it's probably better to keep them in their native UTM projections. Should you wish to ...


1

You are retrieving the coordinates correctly. Google Earth requires all data to be EPSG:4326 - WGS84 Geographic. NOAA converts this data for the KML files however the shapefiles are in an esoteric projection which QGIS identifies as "Unknown datum based upon the Authalic Sphere - Projected" for the shapefiles on that page. The map halfway down the ...


1

You can achieve this using overlay operations. Here's a quick example using some fake data. import geopandas as gpd from shapely.geometry import Polygon # Creating the GeoDataFrame with the grid geometries grid_gdf = gpd.GeoDataFrame(data={'grid_id':[101,102,103,104], 'grid_cat':['W','X','Y','Z'], ...


1

A quick way is to use processing in conjunction with QGIS expressions: using aggregate algorithm, with the trick of using a unique-valued attribute in Group by expression: use the field calc icon to add the desired expressions, this can be done in all fields. output:


1

An acceptable solution for me is to create a virtual layer and then make it permanent/save it. This solution allows me to use the SQL statements from my collection. In a next step I can load them directly from the text files and iterate over a list of queries to automate the process. Furthermore I do not have to touch my original shapefiles. Anyway, I would ...


1

It's often easier to use the GDAL python wrappers. See my attempt below. If you still get errors, check that both raster and shapefile have geographic information (i.e. you can overlap them in QGIS or something like that). from osgeo import gdal from osgeo import ogr def dump_poly(raster_fname, vector_fname, ifeat): # New filename. Assumes input raster ...


1

Use Spatial Join. inner will keep only the grids which overlap polygons [['geometry']] is to select only the polygon geometry column, otherwise you would get all attributes from the polygons in the result. import geopandas as gpd grid = gpd.read_file(r'C:\GIS\data\tempdata\grid.shp') polygons = gpd.read_file(r'C:\GIS\data\tempdata\my.shp') ...


1

This is the public WMS URL for "govMap Israel": https://open.govmap.gov.il/geoserver/opendata/wms? However, there are only few layers which are available from that service, as listed in the screenshot from QGIS software:


Only top voted, non community-wiki answers of a minimum length are eligible