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