texts: OpenLayers 3.5 WFS-T update to GeoServer 2.7 raises StringIndexOutOfBoundsException

texts:
I am using ol.format.WFS.writeTransaction() (OL3.5) to generate WFS-T statements. This works well for insert and delete. The problem is with the update statement.

Any update statement returns following exception from GeoServer:

<ows:Exception exceptionCode="NoApplicableCode">
<ows:ExceptionText>java.lang.StringIndexOutOfBoundsException: String index out of range: -1 String index out of range: -1</ows:ExceptionText>
</ows:Exception>

Edits via QGIS work as expected.

This is an example POST generated by OL3.5

<Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <Update typeName="feature:playa_poly_noattribute" xmlns:feature="http://argeomatica.com">
    <Property>
      <Name>boundedBy</Name>
      <Value>-9693316.43121855,2347365.748407976,-9692555.877885746,2347961.5151853375</Value>
    </Property>
    <Property>
      <Name>geometry</Name>
      <Value>
        <Polygon xmlns="http://www.opengis.net/gml">
          <exterior>
            <LinearRing>
              <posList>-9693072.420357607 2347961.5151853375 -9693316.43121855 2347628.7731022364 -9692822.071552228 2347365.748407976 -9692555.877885746 2347758.700963257 -9692555.877885746 2347758.700963257 -9692737.951605482 2347825.076779921 -9692893.214319186 2347748.042587353 -9693072.420357607 2347961.5151853375</posList>
            </LinearRing>
          </exterior>
        </Polygon>
      </Value>
   </Property>
   <Filter xmlns="http://www.opengis.net/ogc">
     <FeatureId fid="playa_poly_noattribute.1"></FeatureId>
   </Filter>
 </Update>
</Transaction>


texts:
I have table that containing road network. Table structure is

CREATE TABLE final_all_india_road_network
(
  gid serial NOT NULL,
  "name" character varying(254),
  code character varying(254),
  lane character varying(254),
  width integer,
  "type" character varying(254),
  subclass character varying(254),
  surface character varying(254),
  oneway character varying(254),
  bridge integer,
  maxspeed integer,
  the_geom geometry,
  CONSTRAINT final_all_india_road_network_pkey PRIMARY KEY (gid),
  CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
  CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 4030)
)

I want to add check geometry constraints to this table. For that I used following Query

"ALTER TABLE final_all_india_road_network
  ADD CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL);"

But while running this Query following error occured

"ERROR:  check constraint "enforce_geotype_the_geom" is violated by some row"


texts: How to calculate distance matrices for points based on common denominator?

texts:
I am dealing with a vector layer containing ~3400 points. All of them have an individual identifying code (Name) and another one (Network) they share with others of the dataset since they belong a larger network. Thus the attribute table looks somewhat like this:

 Network  Name   AdditionalInformation
 a        Name1  Some
 a        Name2  Information
 a        Name3  Is
 a        Name4  Written
 b        Name5  In
 b        Name6  This
 b        Name7  Field

What needs to be done is calculating the maximum distance between several points belonging to the same network (these are more than 500). However, my knowledge of calculation tools of this kind is limited to QGIS' Distance Matrix tool and I do not see a suitable solution for this entire dataset.

What might be the best way to proceed here?

Is there any way, for instance, to use the Python console in order to calculate distance matrices for each and every network based on the common denominator Network and save all maximum values to a CSV file afterwards?

I split the layer into smaller ones and saved them accordingly. The script, copied from Joseph, is now:

import glob, os, processing

path = "G:\ZIM\Analyse_Foerderprogramme\_ZIM-SOLO und ZIM-KOOP 2014\GIS\Adressen_TOMTOM\Split"
output_path = "G:\ZIM\Analyse_Foerderprogramme\_ZIM-SOLO und ZIM-KOOP 2014\GIS\Adressen_TOMTOM\Split\Ergebnisse"
field_1 = "NAME"
field_2 = "GEOCODE"

for layer in glob.glob(path + "*.shp"):
    processing.runalg("qgis:distancematrix", layer, field_1, layer, field_2, 2, 0, output_path + os.path.basename(layer) + ".csv")

Yet nothing is happening so far after typing it into the console. I am but a beginner in the usage of the Python console so surely I made a mistake. What next?

The earlier mistake was the use of \ instead of the proper /. Now, the script (with more readable paths) looks like:

import glob, os, processing

path = "H:/Split/"
output_path = "H:/Split/Ergebnisse/"
field_1 = "NAME"
field_2 = "GEOCODE"

for layer in glob.glob(path + "*.shp"):
    processing.runalg("qgis:distancematrix", layer, field_1, layer, field_2, 2, 0, output_path + os.path.basename(layer) + ".csv")

This results in the error [Errno 2] No such file or directory: 'H:/Split/Ergebnisse/Adressen_projiziert_Info_NAME__G533.shp.csv' See log for more details

Why does Python look for a CSV file at the very beginning instead of calculating first and writing the results in a file of this sort second?

Joseph's script works. However, the error message remains the same for the file I'm trying to work with. Can anyone tell me what the log message might point to?

2016-09-23T16:22:45 2   Uncaught error while executing algorithm
            Traceback (most recent call last):
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\core\GeoAlgorithm.py", line 203, in execute
                self.processAlgorithm(progress)
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\algs\qgis\PointDistance.py", line 89, in processAlgorithm
                self.writer = outputFile.getTableWriter([])
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\core\outputs.py", line 245, in getTableWriter
                return TableWriter(self.value, self.encoding, fields)
              File "C:/PROGRA~1/QGIS2~1.14/apps/qgis-ltr/./python/plugins\processing\tools\vector.py", line 668, in __init__
                with open(self.fileName, 'wb') as csvFile:
            IOError: [Errno 2] No such file or directory: 'C:/Users/klinge/GIS/Split/Ergebnisse/Ausgangsdaten_Vorhaben__ZKN060503.csv'

2016-09-23T16:22:45 1   There were errors executing the algorithm.

The error message did not occur on another (Linux) computer.


texts: Running multiple iterations of ArcMap model

texts:
I am working with a model that is build in modelbuilder and I have to test with with over 150 different use cases. Is there a way I could upload a table with use cases as it attribute and then have the model run one time on all the 150 cases, one by one. For instance, have it run the first time one case #1 and after is has finished run again at case #2 and so on. I know I have to use iteration but unfortunately it only runs ones.

Below is an image of the model:


texts: Open source space time series database management system

texts: I am aware of the existence of space time DBMSs such as Druid, InfluxDB or OpenTSDB. Is there anything like it for spatially explicit data? A time series extension for a common spatial DBMS (say, PostGIS) could do it; as well as a spatial extension to an existing time series DBMS.


texts: OpenLayers 4 ol.View extent doesn't work horizontally

texts:
I'm attempting to use the ol.View's extent property to ensure the user stays within the bounds of the map. The code below works for the height (vertically) but he user can still move the map outside width wise (horizontally) and see the edges:

var width = window.innerWidth;
  var height = window.innerHeight;
  var mapW = 1789;
  var mapH = 993;

  console.log(width);
  console.log(height);

  var extent = [0, 0, 1789, 993];
  var projection = new ol.proj.Projection({
    units: 'pixels',
    extent: extent
  });

  var map = new ol.Map({
    layers: [
      new ol.layer.Image({
        source: new ol.source.ImageStatic({
          url: 'assets/img/map-bigger.jpg',
          projection: projection,
          imageExtent: extent
        })
      })
    ],
    target: 'map',
    view: new ol.View({
      projection: projection,
      center: ol.extent.getCenter(extent),
      zoom: 3.5,
      maxZoom: 4.5,
      minZoom: 3,
      extent: extent
    })
  });

  map.setView(
    new ol.View({
      projection: projection,
    center: ol.extent.getCenter(extent),
    zoom: 3.5,
    maxZoom: 4.5,
    minZoom: 3,
    extent: map.getView().calculateExtent(map.getSize())
    })
  );

If I resize the window's width to offset this gap it works as expected. However, if I reload the page after the width change the map won't drag far enough to see the full map, as it's cut off.


texts:
I use the tool osmupdate to update the pbf file I get from http://download.gisgraphy.com/openstreetmap/pbf/ (without a pbf suffix by default).

At first merging change files can be finished smoothly but soon it got a error like this:

Killed
osmupdate Error: Merging of changefiles failed:
osmconvert --merge-versions  "osmupdate_temp/temp.d000001399.osc.gz"
"osmupdate_temp/temp.d000001400.osc.gz"
"osmupdate_temp/temp.d000001401.osc.gz"
"osmupdate_temp/temp.d000001402.osc.gz"
"osmupdate_temp/temp.d000001403.osc.gz"
"osmupdate_temp/temp.h000033666.osc.gz"
"osmupdate_temp/temp.h000033667.osc.gz" "osmupdate_temp/temp.8"
-- timestamp=2016-07-16T08:07:02Z --out-o5c >"osmupdate_temp/temp.9"

I have found that the problem lies in shortage of memory since the temp data is much bigger than the original data. Since I only need the updated nodes, how can I get the updated nodes information only?


texts: Distance to inland boundaries rather than coasline

texts:
I have a continent map of Africa with country boundaries. I have converted polygons to lines so that I can calculate the shortest distance of each point (city for instance) to the country border.

However, I want the v_distance function in Grass plug-in (QGIS) to take into account only distance to inland borders, not taking into account the see border.

Do you have an idea of how to perform this ?


texts: I am using QGis 3.2 and on the edit bar I don't have the move feature icon available (even before loading the file). I have seen situation where it is blocked but without showing completely on the bar is a novelty to me. Can anyone give me a hand?


texts: Converting WFS (from ArcGIS Server) to KMZ or SHP

texts:
I'm trying to create a way users can select features from our WFS that's running on ArcGIS Server and have it convert them in KMZ or shapefile. We will be using WebApp builder. I'm pretty comfortable with creating a widget or publishing a geoprocessing tool. I'm good with Python and JavaScript.  I'm just wondering what is the best approach. I looked into FME Server, but that costs a lot of money. I'm wondering if there is a way to do it with GDAL - org2org.
The way we do it now, is the user has to select the features with the select tool in Web AppBuilder. Export it as a GeoJSON and then ship it off to a 3rd party called https://geoconverter.hsr.ch/ . We really do not want to be sending our data there.
Any ideas on how I can do this?
I also tried the clip and ship method Esri has recommend but still it will not convert to KMZ unless you have the iterability extension.  https://enterprise.arcgis.com/en/server/10.3/publish-services/linux/gp-service-example-clip-and-ship.htm


texts: How to read pixel values from an image server layer with ArcObjects?

texts:
I'm having a problem with reading pixel values from the image server layer. The image server layer class implements IRasterLayer, so when I tried to read pixel values using IRaster2 obtained through the IRasterLayer implementation, I could read pixel values. However, these pixel values are sometimes one pixel shifted because there are two rasters in the layer and their pixels are not lined up perfectly.

IRaster2 only accepts one row & column pair, which means that it doesn't take into account multiple raster sources in the image server layer and, if their pixels are not aligned perfectly, just one row & column counting from the upper left corner of the layer cannot read the original pixel values from the data sources.

I tried to access individual rasters stored in the attribute table of the image server layer (Raster column), but it returns null. For example, the following line returns null:

((IAttributeTable)imageServerLayer).AttributeTable.GetRow(id).Value[rasterFieldIndex]

In ArcMap, I can access each row's raster data in the attribute table.

The Identify tool in ArcMap doesn't have this problem and it can read original pixel values from individual raster sources.

I'm sure there should be an easy way of accessing source rasters in the image server layer. Or there should be a recommended method for reading pixel values from the image server layer.


texts: vue.js base wont load in while loading OpenLayers map

texts:
While using OpenLayers library in vue.js the map won't load. It isn't showing any error but only zoom in, zoom out and OSM atrribution is showing.
<template>

</template>


texts: How to split multiple polygon layer into multiple layers by vertices parts?

texts: I have a layer in Arcmap 10 that was brought over from Auto CAD (layer A).  This layer has multiple polygon shapes.  I need to split these shapes from layer A, and place them into their own layers.  When I brought over the layer from Auto CAD, it put all the polygons into one shape.  So when I go to select one in editing, it selects them all.  When I look at editing the vertices I can see each one is a separate part.  Is there a way to split these parts apart to place into their own layers?


texts:
I have two really large datasets, one with about 3 million hexagon features that span the extent of Pennsylvania being used as 'zones', and one with about 4 million features.

I need to run a tabular intersect with these datasets and then perform some formatting and deletion of records. My script works on a subset of the data (at the county level). When I try to run on the full datasets, it runs for about two hours and then produces an empty dataset.  This is understandable because I imagine something is happening with the size of the datasets.

However, I tried to create a loop to select hexagons by county and perform the intersect with those hexagons. When I try this, it will successfully perform the first intersect in the loop and output a table for Adams County, and then fail.  When I ran it in PyScripter, I got the following error:

[Errno 10054] An existing connection was forcibly closed by the remote host

Following the advice here, I uninstalled PyScripter and began just using IDLE as my IDE.  When run from IDE, no error is given, it just ends after completing the first table in the loop. When I copy and paste it into the Python window in ArcGIS, it just crashes ArcMap.  Any ideas?

Here is my code:

import arcpy, os, datetime
from arcpy import env
from arcpy.sa import *
from itertools import groupby
from operator import itemgetter

arcpy.env.overwriteOutput = True

target_features = r'C:\Users\mmoore\Documents\ArcGIS\COA_Species.gdb\PlanningUnit_Hex10acre'
join_features = r'C:\Users\mmoore\Documents\ArcGIS\COA_Species.gdb\SGCN_SDM'
outTable = r'C:\Users\mmoore\Documents\ArcGIS\COA_Species.gdb'
pa_county = r'C:\Users\mmoore\Documents\ArcGIS\COA_Species.gdb\PA_County'

county_lyr = arcpy.MakeFeatureLayer_management(pa_county, "county_lyr")
target_lyr = arcpy.MakeFeatureLayer_management(target_features, "target")

with arcpy.da.SearchCursor(pa_county, "COUNTY_NAM") as cursor:
    county_list = sorted({row[0] for row in cursor})
counties = [x.encode('UTF8') for x in county_list]

for county in counties:
    county_selection = arcpy.SelectLayerByAttribute_management(county_lyr, "NEW_SELECTION", '"COUNTY_NAM" = ' + "'%s'" %county)
    target_selection = arcpy.SelectLayerByLocation_management(target_lyr, "INTERSECT", county_selection, "", "NEW_SELECTION")
    intersect = os.path.join(outTable, county)
    classFeatures = ["ELCODE", "SeasonCode", "OccProb"]
    arcpy.TabulateIntersection_analysis(target_selection, "unique_id", join_features, intersect, classFeatures)

    case_fields = ["unique_id", "ELCODE", "SeasonCode"]
    max_field = "PERCENTAGE"
    sql_orderby = "ORDER BY {}, {} DESC".format(",".join(case_fields), max_field)

    with arcpy.da.UpdateCursor(intersect, "*", sql_clause=(None, sql_orderby)) as cursor:
        case_func = itemgetter(*(cursor.fields.index(fld) for fld in case_fields))
        for key, group in groupby(cursor, case_func):
            next(group)
            for extra in group:
                cursor.deleteRow()


texts: GeoServer: Receive multiple tiles from layer

texts:
Using GeoServer you can do a WMS request to receive a single PNG representation of a certain layer using something like this.

http://localhost:8080/geoserver/wms?service=WMS&version=1.1.0&request=GetMap&layers=spearfish&styles=&bbox=589425.9342365642,4913959.224611808,609518.6719560538,4928082.949945881&width=512&height=359&srs=EPSG:26713&format=image%2Fpng

Is it possible to alter this request to directly request different "zoom levels" or "tiles" of the layer? I assume you need to do a Web Map Tile Service (WMTS) request but I can't seem to find any documentation on how to do this. Any guidance would be helpful.


texts: Using list of PIN14 numbers and addresses to replace addresses in ArcMap attribute table?

texts:
In my attribute table I have a PINs field and an Address field (along with other zoning information) that I want to replace with correct addresses.

Is there a way I can have ArcMap go through and replace all the addresses with the new ones so long as they match an existing pin and if there is no match then have it remove the current address? I assume there is some simple field calculator expression for this.


texts: Is there an easy way of get related table attributes to appear all on the first page of a pop up in an ArcGIS online web map?


texts: Classification of multiple values at the same point

texts: If I have a point dataset with multiple values for the same geographic location (e.g. residuals from a regression where multiple individuals live at the same location) which I classify using quantile classification in QGIS, how does QGIS choose which quantile (in practice which color) to show on the map? For example, if there are three values for a single location with values falling within e.g. the first, second and fifth quantile, how does QGIS choose what to represent on the map?


texts: In QGIS, how to input & select multiple values from GUI? (without repetitive OR, OR, OR)

texts:
In QGIS, most likely from the SELECT FEATURES BY VALUE...(F3) or SELECT FEATURES BY EXPRESSION...(CTRL-F3), what is the simplest most flexible way to put in multiple values and have QGIS show that selection on the canvas?
I have list of comma separated values which are the values of features I want to select from a spatial dataset. 1, 2, 3, 4, 5. Because it might be a series of 50 values, I want to quickly paste that into an expression, without adding single quotes or the like.
I'm thinking array_to_string is part of the story, but not seeing any example of how that goes together. Unnest doesn't seem to be a supported function. My test below.
Or possibly in the DBManager as SQL, however I think my colleagues will not adopt that as readily.


texts: Identify circuits/loops in GPS track data

texts:
I have a GPS time series with a set of lat/long coordinates per second, how would you go about identifying circuits/loops in that GPS track?

To the layman, it seems like something that should be possible, but I suspect it could just as well be a tremendously complicated problem.


texts:
Is there a simple way to save a QGIS project and all associated files to a zip folder so that another person could unzip the folder and open the project?

Or is there a better way to share a project with another person?


texts: Saving individual figures in Print Composer in QGIS 2

texts:
Can anyone tell me how to save individual figures in QGIS print composer that I can later go back and edit?

Do I have to create a new print composer for each figure?

I am using QGIS 2.18.9.


texts: Is it possible to edit Delimited Text Layers in QGIS?

texts:
GIS beginner here...

I have a tab-delimited text file that I'm importing into QGIS just fine. But I try to edit the attribute table, the little pencil-looking button is grayed out, and it's not letting me edit it via the QuickMultiAttributeEdit plugin either.

I've googled around a lot, but can't for the life of me figure out how to make the attribute table editable. Any help?

Thanks all!


texts: How to get X Y coordinates and cell value of each pixel in a raster using Python?

texts:
I'm really new to Python and I would like to know whether there is a quick method to get cell values of a raster pixel by pixel and the coordinates (map X Y coordinate of centre of each pixel) using Python in ArcGIS 10?

To describe this further, I need to get the map X, map Y and cell value of the first pixel and assign those three values in to three variables and repeat this step for rest of the other pixels (loop through the entire raster).

I think I need to describe my question more. Problem is, I need to get the X Y location of a pixel of the first raster and get cell values of several other rasters corresponding to that X Y location. This process should be loop through every pixel of the first raster without creating any intermediate point shapefile as it is going to be really really time consuming since I have to handle a raster with nearly 8 billion pixels. Also, I need to do this using Python in ArcGIS 10.

@JamesS: Thank you very much for your suggestion. Yes this would work for one raster but I need to collect the cell values for several other rasters as well. The problem is, after getting the X and Y coordinate of the first pixel of the first raster, I need to get the cell value of the second raster corresponding to that X, Y location of the first raster, then third raster and so on. So, I think when looping through the first raster, getting X and Y location of a pixel and getting cell values of the other raster corresponding to that location should be done simultaneously but I'm not sure. This can be done by converting the first raster in to a point shapefile and performing Extract multivalues to point function in ArcGIS 10 but I'm unable to follow that method because I don’t want to create a shapefile as it is going to be really slow and I HAVE TO FIND A SOLUTION USING PYTHON, not with any existing tool in ArcGIS.

@hmfly: Thanks, Yes this method (RastertoNumpyarray) will work if I can get the coordinate of a known row and column value of the array.

@whuber: I don’t want to perform any calculations, all I need to do is write X Y coordinates and cell values in to a text file and that’s all


texts: RasterIO saves wrong coordinate system when writing to a new file

texts:
I am trying to open, perform a calculation then save to a new tif file using RasterIO without changing the coordinate system.

My input tif is in GDA 1994 MGA Zone 54 and has a wkt string of

 PROJCS["GDA_1994_MGA_Zone_54",GEOGCS["GCS_GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

RasterIO, for the same tif,  reads and writes the coordinate system as  UTM Zone 54, Southern Hemisphere and which has a wkt string of

PROJCS["UTM Zone 54, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]

Both these coordinate system have the same Proj4 string. It seems that RasterIO creates its crs.wkt string from the proj4 string, not from the one provided by the input tif.

Is there a way using RasterIO to read the wkt string from the tif and use it when saving to a new file?

Here's my code so far:

smin=0
smax=255
with rasterio.open('input_94mga54.tif') as src:
    band = src.read(1,masked=True)
    rescaled = (band - band.min()) * ( smax - smin) / (band.max() - band.min()) + smin
    with rasterio.open('rescaled_94mga54.tif', 'w', **src.meta) as dst:
        dst.write_band(1, rescaled)


texts: QGIS (merging, dissolve, union) lines of same name

texts:
I have street data where roads are segmented for no particular reason. I want to merge, dissolve or union whatever seems fit to make a single flat line. I H
have tried above processes in QGIS but every time I start edit in QGIS and start deleting the vertexes of that line the vertex jumps around as if that line is still segmented. I have attached one image where you can see how a dissolved line when edited its vertex are not deleted in straight path rather its jumps around.
So the questions is how can I make a single, straight, flat line out of segmented based on road names.


texts: Is it possible to save qgs file as older version?

texts: I just got installed the version 2.6 of QGIS (OSGeo4W) on Windows 7. I am exploring the new capabilities and the compatibility of our tools with this new version, however, i have to keep on working on the old gqs files, because we have to send them to other users who do not have the newest versions, so, I must to be able to save my qgs files as 2.2 Valmiera version, is this possible?


texts:
I have few files created with ArcGIS in .mxd format.

Now I dont have ArcGIS and want to work with QGIS along with these .mxd file formats.

How can I solve this problem (along with each step to be taken)?


texts: How to delete geohash rectangles from Leaflet?

texts:
I have been inspired by this code to create a geohash map for my Neo4j geohash graph on Leaflet:

Unlike the (lat/lon) markers, the geohash string is represented as a rectangle in the map.

For a static search, I return all the geohash nodes with an Ajax query in my Leaflet Map.

For a dynamic search I added a search  form, as you can see above the map, to query the Neo4j database by time (timehash here for me).

The request to the database is working fine, and the new geohash rectangles are plotted in the map, however the previous geohash rectangles remains.

How do I delete the previous geohash rectangle from the map?


texts:
I need to compute the centroid position of an irregular polygon with lat and lon coordinates in Java.
Which formula should I use to correctly compute a centroid?


texts: Setting QGIS simple line color depending on fields value in python plugin in QGIS

texts:
I am new to QGIS and using QGIS 3.6. I am trying to add a geojson file to QGIS and then export it to a png.

I have a color filed in geojson, and can apply it with GUI.
Now I want to do it using python plugin.
In the GUI, I can easily go to symbology and choose simple line, set fields as color.

I tried some code in plugin

layer = iface.addVectorLayer("/Users/usr/Downloads/test8.json", "my layer", "ogr")
layer = iface.activeLayer()
single_symbol_renderer = layer.renderer()

symbol = single_symbol_renderer.symbol()
symbol.setColor(QColor.fromRgb(0, 225, 0))
symbol.setWidth(0.75)

layer.triggerRepaint()

The code above can change the color to green. But I could not find any way to set the color using the color field.


texts: Harmonic modeling of multiple bands in Earth Engine

texts:
I am trying to generate harmonic models for multiple spectral indices derived from the same image collection time-series in Google Earth Engine (say EVI and NDVI).
The implementations I have been able to find of harmonic time-series analysis (for example here, here, and here, all of which seem to be adaptations of Nicholas Clinton's 2016 example) all model a single variable using multiple, sequential functions. I am struggling to figure out how this could be adapted into a single function, such that I could just map() over an image collection and produce a harmonic model for each band of interest. My first thought was just to provide the band names as a list but that doesn't work because the linearRegression reducer doesn't accept lists. And the temporal segmentation algorithms seem like they do what I need but then also do a bunch of other stuff that I don't need...namely temporal segmentation.
This seems like it has to be pretty simple, but I'm stuck on it.
Example code (not my own! This is from the good people at SERVIR Mekong) illustrating what I mean about sequential functions is copied below.
// The dependent variable we are modeling.
var dependent = 'EVI';

// The number of cycles per year to model.
var harmonics = 3;

// Make a list of harmonic frequencies to model.
// These also serve as band name suffixes.
var harmonicFrequencies = ee.List.sequence(1, harmonics);

// Function to get a sequence of band names for harmonic terms.
var getNames = function(base, list) {
  return ee.List(list).map(function(i) {
    return ee.String(base).cat(ee.Number(i).int());
  });
};

// Construct lists of names for the harmonic terms.
var cosNames = getNames('cos_', harmonicFrequencies);
var sinNames = getNames('sin_', harmonicFrequencies);

// Function to add a constant band.
var addConstant = function(image) {
  return image.addBands(ee.Image(1));
};

// add a constant as a band
var harmonics = EVI.map(addConstant);

// Function to add a time band.
var addTime = function(image) {
  // Compute time in fractional years since the epoch.
  var date = ee.Date(image.get('system:time_start'));
  var years = date.difference(ee.Date('1970-01-01'), 'year');
  var timeRadians = ee.Image(years.multiply(2 * Math.PI));
  return image.addBands(timeRadians.rename('t').float());
};

harmonics = harmonics.map(addTime);

// Function to compute the specified number of harmonics
// and add them as bands.  Assumes the time band is present.
var addHarmonics = function(freqs) {
  return function(image) {
    // Make an image of frequencies.
    var frequencies = ee.Image.constant(freqs);
    // This band should represent time in radians.
    var time = ee.Image(image).select('t');
    // Get the cosine terms.
    var cosines = time.multiply(frequencies).cos()
      .rename(cosNames);
    // Get the sin terms.
    var sines = time.multiply(frequencies).sin()
      .rename(sinNames);
    return image.addBands(cosines).addBands(sines);
  };
};

harmonics = harmonics.map(addHarmonics(harmonicFrequencies));

print(harmonics)

// Independent variables.
var independents = ee.List(['constant', 't'])
  .cat(cosNames).cat(sinNames);

// The output of the regression reduction is a 4x1 array image.
var harmonicTrend = harmonics
  .select(independents.add(dependent))
  .reduce(ee.Reducer.linearRegression(independents.length(), 1));

// Turn the array image into a multi-band image of coefficients.
var harmonicTrendCoefficients = harmonicTrend.select('coefficients')
  .arrayProject([0])
  .arrayFlatten([independents]);

// Compute fitted values.
var fittedHarmonic = harmonics.map(function(image) {
  return image.addBands(
    image.select(independents)
      .multiply(harmonicTrendCoefficients)
      .reduce('sum')
      .rename('fitted'));
});


texts: Does anyone know the difference in disable 'Cache Feature' and 'Only request features overlapping the current view extend'. The result seems to be the same. I would expect the disabling of 'Cache Feature' to dynamically load WFS features when panning or zooming the map.


texts: Batch processing - Translate - JPEG compression parameters

texts:
How can I configure the batch processor to enable JPEG compression?
Just like in the "Run as single process..."
I guess I need to add instructions "Additional command-line parameters".


texts: Spatial downscaling of MODIS LST grids using a covariate (e.g. NDVI): available software?

texts:
I have a MODIS MOD11A2 (LST) grid, which has a spatial resolution of 1000m. For the purpose of Surface Energy Balance modelling, I'm trying to improve the resolution of this data to 250m, using a covariate like the NDVI. Lots of literature is available on this topic, with various titles and subjects like "spatial downscaling', 'downscaling co-kriging', etc.

I've been looking for software that is capable of performing this operation for me. Since I'm a R user, I was happy to find the dissever package (https://github.com/pierreroudier/dissever). However, this package seems like it is no longer maintained, and some errors prohibid me from using the package.

Does anyone know a working sofware package, preferably in R or Python, that is capable of performing such spatial downscaling?


texts: UK Ordnance Survey Code Point. Projection change (OSGB 1936 to WSG84)

texts:
I am trying to overlay CodePoint with Polygons onto leaflet JS. My current problem is converting the .shp file into the correct projection and a kml.

If I open the shp file (using qgis) with the projection of "OSGB 1636" which (from what I understand is the correct projection), when I "Save as" and select the new projection of "WSG84". This is the result:

  <Placemark>
    <Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
  <MultiGeometry><Polygon><outerBoundaryIs><LinearRing><coordinates>inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf</coordinates></LinearRing></outerBoundaryIs></Polygon><Polygon><outerBoundaryIs><LinearRing><coordinates>inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf inf,inf</coordinates></LinearRing></outerBoundaryIs></Polygon></MultiGeometry>

Using ogr2ogr command:

ogr2ogr -t_srs 'EPSG:4326' -s_srs 'EPSG:4277' -f 'KML' 84.kml SO.shp

I get this output:

Failed to reproject feature 0 (geometry probably out of source or destination SRS).
ERROR 1: Terminating translation prematurely after failed
translation of layer SO (use -skipfailures to skip errors)

Is there anything I'm doing incorrectly - to get this data to the correct projection?


texts: How can I change the frame in the print composer in QGIS?

texts:
I made a map based on a raster from Natural Earth web site with some vector layers with the Qgis software.
When I compose this map I cannot get the right size: there are too big or to small as you can see in the photos. I tried to change all the setting in the composer, but does not work.

Can anywhere, please, can explain how to crop the frame in the composer so I can get rid of the white border?


texts: Referencing layer in table of contents using PyQGIS?

texts:
I'm trying to create a plugin for QGIS using python, But need some help. I want to 1.Zoom to a layer and 2. make the layer visible(switch off or on).

But I have managed to zoom to an active layer by setting the active layer in the map canvas by using.

import qgis
canvas = qgis.utils.iface.mapCanvas()
acl = canvas.layer(1)
qgis.utils.iface.setActiveLayer(acl)
> True
qgis.utils.iface.zoomToActiveLayer()

But this only work as long as the layer is viable(obviously since I use canvas). So I need another way to "point" to a layer and then zoom to and switch off.


texts: Zooming to view more/less details of wms layer in Leaflet

texts:
I have a layer in geoserver that contains data of different levels: Country boundary, Province Boundary and District boundary in a single file(distinguishable from attribute fields). I have loaded the geoserver layer in Leaflet as tiled wms layer. I want to show only country boundry at first and when the user zooms in to an extent (say zoom level 7), then show the provincial bounds and so on.Is it possible to do so with only one layer considering attributes? (I don't want to load multiple layers).
This is what I've done to load wms layer:

var mywms = L.tileLayer.wms("http://localhost:8080/geoserver/topp/wms", {
    layers: 'topp:boundaries',
    format: 'image/png',
    transparent: true,
});

mywms.addTo(map);


texts: Saving layer as temporary file usig PyQGIS?

texts:
I am writing a script to automate some processes which involved producing some layers like this:

Input -> Do sth on Input -> Layer A -> Do sth on Layer A-> Layer B -> Do sth on Layer B -> output

My approach is to save Layer A, B and the output as shapefile

However, layer A B are actually useless after the process is finished. I don't want to save them as shp but as temporary layer like below:

But I want to do it via PyQGIS


texts: How to make an emf into a character marker symbol north arrow?

texts: I have an emf that I would like to make into a North Arrow. How do I get the emf to be a Character Marker Symbol in ArcGIS 10?


texts:
An institution wants to have layer/document uploads done only by specific users.
The "normal" user shall be restricted from uploading a layer/document.

Can this be achieved in the recent version of geoNode?
If yes, how?


texts: Create a list and append objects created from function arcpy.Exists through a loop

texts:
I used a SearchCursor object to create a list
 containing file paths and saved the list to a variable called FilePathList.

Next, I iterated through elements within the list using a for loop. Within the for loop, I nested an if statement that contains the arcpy.Exists function. I am inserting elements from the FilePathList into the arcpy.Exists function. When list elements are inserted into the function, if the element exists a data object string “1” is created. In addition, if the element does not exist a data object string "0" is created. I read that the arcpy.Exists() function returns a Boolean indicating whether the element exists.

I want to use the arcpy.Exists function to create string objects "1" representing exists and "0" representing does not exist. In the code below, I am able to print the string objects one by one, in the correct order, to the console. However, I would like to append the strings to an empty list in the order that they are printed. The order is important because I will eventually insert the created list into a new field within the dbf file. The field needs to be in the correct order because it corresponds and links to another field.

import arcpy

arcpy.env.workspace = r"C:\TestFolder"

#dbf table used for processing
dbf = "C:\\TestFolder\\ProcessControlTable.dbf"
#Column within dbf table containing file paths
myColumn = "FI_PATH"
#Temp list created to hold appended items from arcpy.Exists Function
Temp = []

#SearchCursorList created by using SearchCursorObject to create List
FilePathList =[row[0] for row in arcpy.da.SearchCursor(dbf,myColumn)]
print(FilePathList)
for elements in FilePathList:
    if arcpy.Exists(elements):
        #If element in SearchCursorList Exists;
        #Define a variable called Exists and string object
        Exists = "1"
        Temp.append(Exists[0])
        print(Exists)

    else:
        DoesNotExist = "0"
        Temp.append(Exists[0])
        print(DoesNotExist)

print (Temp)

The output in the Python shell first prints the list containing file paths.
Next, strings created from looping with the arcpy.Exists function are printed one by one in the correct order. Lastly, the Temp list is printed however all items in the list are the string object '1'. I want Temp list to contain strings created from the arcpy.Exists function looping in the correct order. The Temp list should print ['1','1','0','1','1','1','1','1','1','0']

I'm not sure what code will pull objects out of the loop and put it into a new list. I can't think of the logic with my limited Python and ArcPy knowledge.


texts:
For some reason I am unable to import shape files in the form of point data to REDCAP.  Do these data points need projected coordinate systems? I am aware that the sample size is 500 when double clicking the REDCAP file.

I am also encountering errors when using the command line to launch the application.

Java is not recognized as an internal or external command when trying to execute the code included in the manual to launch the application.  There are apparently memory unless launched from Java.  This enables working with larger data sets.

Neither polygons or points are opening, the data have been tested in a GCS and PCS.
The sample data opens without flaw.

http://www.spatialdatamining.org/software/redcap

Any suggestions would be much appreciated.


texts: How can I remove small polygons from a raster data set? I know that they aren't technically polygons, but I convert the rasters to polygons. The problem is that I want to focus on a specific polygon that is large and there are many small polygons surrounding it. I want the small polygons out. Can anyone help me with this?


texts: PyQGIS3 Layout : item position relative to another item size

texts:
I created an atlas layout with a map, an attribute table, and a label.
I would like to make the y position of the label move according to the size of the attribute table on each atlas page.
The code below allows to display these three elements on the layout, but the label position is determined only according to the size of the attribute table on the last page.
from collections import *

project = QgsProject.instance()
manager = project.layoutManager()
layoutName = "TEST"

layouts_list = manager.printLayouts()
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)

layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(layoutName)

myAtlas=layout.atlas()
myAtlas.setFilterFeatures(True)

manager.addLayout(layout)

layout.atlas().setEnabled(True)

pc = layout.pageCollection()
pc.page(0).setPageSize('Letter', QgsLayoutItemPage.Orientation.Landscape)
cover = project.mapLayersByName('client_Atlas')[0] #os.path.basename(path) +
myAtlas.setCoverageLayer(cover)
myAtlas.setHideCoverage(True)
myAtlas.setSortFeatures(True)
myAtlas.setSortExpression('diagramme')
myAtlas.setFilterExpression('diagramme')

""" Map"""
map = QgsLayoutItemMap(layout)
map.setId('map1')
map.setRect(20, 20, 20, 20)#Set Map Extent
map.setAtlasDriven(True)
map.setAtlasMargin(0.2)
rectangle = QgsRectangle(-244516.858, 169608.638, -243208.118, 170778.898)
map.setExtent(rectangle)
layout.addLayoutItem(map)
map.attemptMove(QgsLayoutPoint(2.7, 30.830, QgsUnitTypes.LayoutMillimeters))
map.attemptResize(QgsLayoutSize(211, 171.8, QgsUnitTypes.LayoutMillimeters))

""" Attribute table """
pdf = project.mapLayersByName('plan_de_ferme')[0]
pdf_table = QgsLayoutItemAttributeTable.create(layout)
pdf_table.setVectorLayer(pdf)
pdf_fields = ['champ','sup','supepnd']
pdf_table.setDisplayedFields(pdf_fields, True)

pdf_table.setDisplayOnlyVisibleFeatures(True)
pdf_table.setFilterToAtlasFeature(True)
pdf_table.setMap(map)
layout.addMultiFrame(pdf_table)
frame2 = QgsLayoutFrame(layout, pdf_table)
frame2.attemptResize(QgsLayoutSize(60.6, 153.3),True)
frame2.attemptMove(QgsLayoutPoint(213, 30.73, QgsUnitTypes.LayoutMillimeters))
pdf_table.addFrame(frame2)

""" Total text label """
title = QgsLayoutItemLabel(layout)
page_field = cover.fields().indexFromName('sup_sum')
for f in cover.getFeatures(): # get sum values for each page
    Tsupha = "[%attribute( get_feature('client_Atlas' ,'diagramme', attribute('diagramme')),'sup_sum')%]"
    Tsupepndha = "[%attribute(get_feature('client_Atlas','diagramme',attribute('diagramme')),'supepnd_su')%]"
    Tsupac = "[%attribute( get_feature('client_Atlas' ,'diagramme', attribute('diagramme')),'supac')%]"
title.setText('Total : ' + Tsupha + '   ' + Tsupepndha + '   ' + Tsupac)
myFont=QtGui.QFont("Ms Shell Dlg 2", 10)
title.setFont(myFont)
title.setHAlign(Qt.AlignLeft)
title.setVAlign(Qt.AlignTop)
title.setMinimumSize(QgsLayoutSize(110.39, 8.41, QgsUnitTypes.LayoutMillimeters))

field = 'diagramme' # attribute that contains categories
c = Counter([f[field] for f in pdf.getFeatures()]) #get the number of features for each category/page

cover = project.mapLayersByName('client_Atlas')[0]
features = cover.getFeatures()
length = 0
for feat in features : #count the number of features in the cover layer
    length += 1
for i in range (0,length):
    nbrow = c.get(i+1) + 1 # number of row + 1 for the headers
    totalH = nbrow*6 # multiplied by the heigth of rows in layout's attribute table
    Ytab = 30.73 + totalH + 5 # y position + heigth of attribute table + 5 mm
    title.attemptMove(QgsLayoutPoint(213, Ytab, QgsUnitTypes.LayoutMillimeters))
    layout.addLayoutItem(title)

myAtlas.beginRender()
for page in range(myAtlas.count()):
    myAtlas.next()

myAtlas.endRender()


texts: Writing layer to Geopackage using PyQGIS

texts:
I have a PyQGIS script which carries out a number of operations on layers and writes the results to Shape files. However I understand that GeoPackages have a number of advantages (not having multiple files being an obvious one) and would like to change to use them. However I am having trouble writing to them.
Using QgsVectorFileWriter.writeAsVectorFormat I was able to write to a GeoPackage. Once. But on calling the script again I get an error "a file system object (file name} already exists". Fair enough: I added option QgsVectorFileWriter.CreateOrOverwriteFile. But this has no effect - I get the same error message. If I use option QgsVectorFileWriter.CreateOrOverwriteLayer then I can overwrite an existing file, but cannot create the file in the first place: it produces an error 'Opening of data source in update mode failed'.
I have tried both writeAsVectorFormat and writeAsVectorFormatV2 but can get neither to work.
This is with QGIS 3.10.
Why does CreateOrOverwriteFile not cause the file to be overwritten? How does one code this so that it works both on the first and the subsequent times that the script is used?
transform_context = QgsProject.instance().transformContext()
save_options = QgsVectorFileWriter.SaveVectorOptions()
#save_options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteFile
save_options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
#err = QgsVectorFileWriter.writeAsVectorFormat(layer, gpkgFilePath, save_options)
err = QgsVectorFileWriter.writeAsVectorFormatV2(layer, gpkgFilePath, transform_context, save_options)

With regards to the first part of my question: it turns out to be a combination of a timing issue and an erroneous error message. Apparently when you remove a Shape layer QGIS closes the file immediately, but when you remove a Geopackage layer it only does so later (even if it is the only layer in the file). If one then tries to overwrite the file it should produce a "file is open" message, but instead says "file already exists" (which is true but irrelevant).
As far as getting it to work is concerned, it works if I do this:
transform_context = QgsProject.instance().transformContext()
save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
err = QgsVectorFileWriter.writeAsVectorFormatV2(layer, gpkgFilePath, transform_context, save_options)
if err:
    if err[0] == QgsVectorFileWriter.ErrCreateDataSource:
        show('-- does not exist; creating file: ' + gpkgFilePath)
        save_options = QgsVectorFileWriter.SaveVectorOptions()
        err = QgsVectorFileWriter.writeAsVectorFormatV2(layer, gpkgFilePath, transform_context, save_options)
        if err and err[0]!=0:
            show("Error creating file with layer " + layerName + " " + str(err))
            return
    elif err[0]!=0:
        show("Error writing layer " + layerName + " " + str(err))
        return

Though that is awfully clunky - it's an answer, but not a terribly good one. Is there no way of doing this with a single write?


texts: Graphic (Dynamic) Table missing from PDF mapbook export when arcpy script tool used

texts:
When exporting a mapbook (data driven pages) using a Python script tool, the Graphic Tables are missing all rows on the exported PDF pages (but the table header is there).  When exporting the mapbook by running the script independently (using same parameters), the Graphic Table is populated as expected.
The empty Graphic Tables are the only differences between the two methods of running the script, otherwise the PDF mapbook is created as expected.
I am using ArcMap 10.5 and Python 2.7.
I have tried checking out the Production Mapping extension, arcpy.CheckOutExtension("Foundation") (which is required for Graphic Tables) and have tried applying the extension check out before creating the mxd object, after creating the mxd object, and even tried checking out the extension within the ArcPy tool validator.
I've also tried setting higher resolutions and adding ddp.refresh() and arcpy.RefreshActiveView() to see if that helped.  It didn't seem to.
When using the same parameters as the script, the script run independently produces a mapbook PDF export exactly as desired, and running from the tool produces the same, with the exception that the Graphic Tables are empty.
It appears that the MXD is not recognizing the Production Mapping extension when the tool runs the script, but is recognizing that extension when the script runs independently.
My code is as follows:
    mxd = arcpy.mapping.MapDocument(exportmxd)
    pdfPath = os.path.join(outputFolder, "{0}_{1}_{2}_{3}.pdf".format(installation, mapName, selection, todayDateString))
    if os.path.exists(pdfPath):
        os.remove(pdfPath)
    finalPDF = arcpy.mapping.PDFDocumentCreate(pdfPath)
    ddp = mxd.dataDrivenPages
    indexLayer = ddp.indexLayer
    expression = "FEATURE_NAME = '{0}'".format(selection)
    arcpy.SelectLayerByAttribute_management(indexLayer, "NEW_SELECTION", expression)
    ddp.refresh()
    arcpy.RefreshActiveView()
    arcpy.CheckOutExtension("Foundation")
    ddp.exportToPDF(out_pdf = pdfPath, page_range_type = "SELECTED", multiple_files = "PDF_SINGLE_FILE", resolution = 300, image_quality = "BETTER")
    arcpy.CheckInExtension("Foundation")


texts: Determining the first and last point that a polygon intersects

texts:
I have this...

And I want to join those attributes to the polygon automatically:

Because I'll have to replicate this on tones of those polygons.

If this works I would like to know if it still works for non-intersections.


texts: PostgreSQL: insert from multiple select?

texts:
Is it possible in PostgreSQL to insert field values into a table using multiple selects from various other tables? Something like the following SQL (I have over simplified my original problem). I don't want to perform two individual insert operations.

insert into table_1 (name, id)
(select name, id from table_2 limit 1),
(select name, id from table_3 limit 1);


texts: Typing path length using shortest path dijkstra()

texts:
I am now trying to get the length of a shortest path using dijkstra(). I need it to find all the points in my sample, which are located more than 25 km away along roads from district centres. I will loop it, but before I do this I need to get the path's length.

When I enter the following code, it only shows me the shortest path highlighted in red, but not its length. How can I modify it in order to get length in km? (In case it is important, my crs is North_Pole_Lambert_Azimuthal_Equal_Area)

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(2327292.145, -2732901.348)
pStop = QgsPoint(2328327.741, -2732822.68)

tiedPoints = director.makeGraph(builder, [pStart, pStop])
graph = builder.graph()

tStart = tiedPoints[0]
tStop = tiedPoints[1]

idStart = graph.findVertex(tStart)
idStop = graph.findVertex(tStop)
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
if tree[idStop] == -1:
  print ("Path not found")
else:
  p = []
  curPos = idStop
  while curPos != idStart:
    p.append(graph.vertex(graph.arc(tree[curPos]).inVertex()).point())
    curPos = graph.arc(tree[curPos]).outVertex();
    p.append(tStart)
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.red)
  for pnt in p:
    rb.addPoint(pnt)


texts: Setting up ODBC Connection with Access database in MapBasic?

texts:
Good day! i'm trying to setup an ODBC connection in my Access database through mapbasic code, but i don't how should i start this code.

I could make it manually and able to have a link between my tables to work on.

I would like to create a tools that will open and link my tables in access and update it time to time when every time i change something on my mapinfo table.


texts: Calculating a Unique ID on a Subgroup of Polygons

texts:
I would like to assign a unique ID to polygons in a dataset.  The catch to my issue is that I would like polygons that fall within a polygon to have the same unique ID.  My example is related to buildings.  I have many buildings in a dataset.  Each building should have a unique ID, however, buildings that are comprised of multiple polygons (see graphic) should have the same unique ID.  I know of a way to do this through a series of dissolves, and spatial joins, but I would like to do this in a calculation.  My preference is to do this in ArcMap, but QGIS would also be acceptable.

Any thoughts or help would be appreciated.

Thanks...

EDIT:

Here is a graphic of what I am hoping to resolve. All touching polygons have the same ID. This is unique in the entire dataset, not with the building that is shown.


texts:
I have NetCDF files with some datasets with data and two datasets containing lons and lats. I would like to open data layers in QGIS in coordinates that are written in NetCDF's lon/lat datasets. How could I do that?

What I am trying to do is transform NetCDF data to another coordinate system in another grid, then save it as TIFF. I wanted to use QGIS first to see what steps should I make.


texts: gdal Import error: DLL Load Failed when executing gdal_calc.py

texts:
I have the recent version of QGIs/GDAL installed.
If I try to execute gdal_calc.py, I get the following error message:
`Traceback (most recent call last):
File "C:\OSGeo4W64\apps\Python37\lib\site-packages\osgeo_init_.py", line 18, in swig_import_helper
return importlib.import_module(mname)
File "C:\OSGeo4W64\apps\Python37\lib\importlib_init_.py", line 127, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
File "", line 1006, in _gcd_import
File "", line 983, in _find_and_load
File "", line 967, in _find_and_load_unlocked
File "", line 670, in _load_unlocked
File "", line 583, in module_from_spec
File "", line 1043, in create_module
File "", line 219, in _call_with_frames_removed
ImportError: DLL load failed: Das angegebene Modul wurde nicht gefunden.
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "C:\Program Files\GDAL\gdal_calc.py", line 5, in
from osgeo.utils.gdal_calc import *  # noqa
File "C:\OSGeo4W64\apps\Python37\lib\site-packages\osgeo_init_.py", line 32, in
gdal = swig_import_helper()
File "C:\OSGeo4W64\apps\Python37\lib\site-packages\osgeo_init.py", line 31, in swig_import_helper
return importlib.import_module('gdal')
File "C:\OSGeo4W64\apps\Python37\lib\importlib_init.py", line 127, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
ImportError: DLL load failed: Das angegebene Modul wurde nicht gefunden.
E:\UniversitaetHD\WS 2020 - 2021\Seminar GIS\Gruppenarbeit\Dateien>PAUSE`
Seems like those problems to me:
ImportError: DLL load failed with error code when I try to import GDAL from the Python console
gdal ImportError in python on Windows
Unfortunately, none of these solutions work for me.
My system variables are already right configured:

While executing from osgeo import gdal in python, the same error message appears.
I installed the latest GDAL versions from here: https://www.gisinternals.com/release.php


texts:
I'm using Leaflet to visualize a large number of points on a map. The points don't need to be very fancy, I just need to be able to set a radius (doesn't matter if it's pixels or meters, I just need relative size differences), and a color on them.

So I'm thinking that L.Circle and L.CircleMarker are overkill for this purpose. With 7000 points, the map lags quite a lot. Are there any markers that are more light-weight than L.Circle/L.CircleMarker? Some kind of markers that doesn't need complicated radius calculations or a lot of different drawing settings?

If not, any other way I can get better performance with that many markers? I tried with a conditionalLayer, but it's too slow to update, since it has to remove and re-add everything for every zoom or pan (which takes 2-3 seconds every time). It is important that I am able to show all the markers (up to 10 000) when zoomed out.

(Note: I have already activated the L_PREFER_CANVAS setting, which did make a difference, but it's still laggy.)


texts:
My file is

Driver: GTiff/GeoTIFF
Files: 0001.TIF
Size is 8512, 8622
Coordinate System is:
PROJCS["WGS 84 / UTM zone 31N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32631"]]
Origin = (837903.261500865686685,5877808.477226150222123)
Pixel Size = (30.079329235631686,-30.079329235631686)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  837903.262, 5877808.477) (  8d 1'48.79"E, 52d56'35.71"N)
Lower Left  (  837903.262, 5618464.501) (  7d46'39.67"E, 50d37'11.83"N)
Upper Right ( 1093938.512, 5877808.477) ( 11d48'11.42"E, 52d43'20.22"N)
Lower Right ( 1093938.512, 5618464.501) ( 11d21'56.03"E, 50d24'59.35"N)
Center      (  965920.887, 5748136.489) (  9d44'38.52"E, 51d41'24.52"N)
Band 1 Block=8512x1 Type=UInt16, ColorInterp=Gray

After reprojecting it to EPSG:32632 by

gdalwarp -t_srs epsg:32632  0001.TIF 0001_.TIF

I got

Driver: GTiff/GeoTIFF
Files: 0001_.TIF
Size is 9192, 9293
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (413533.214507141499780,5866398.182513756677508)
Pixel Size = (29.999953551392952,-29.999953551392952)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  413533.215, 5866398.183) (  7d42'47.95"E, 52d56'24.88"N)
Lower Left  (  413533.215, 5587608.614) (  7d46'56.92"E, 50d26' 2.84"N)
Upper Right (  689292.788, 5866398.183) ( 11d48'55.18"E, 52d54'49.85"N)
Lower Right (  689292.788, 5587608.614) ( 11d39'51.06"E, 50d24'35.94"N)
Center      (  551413.001, 5727003.398) (  9d44'37.73"E, 51d41'29.87"N)
Band 1 Block=9192x1 Type=UInt16, ColorInterp=Gray

I calculated myself (using Java GeoTools) that

  413533.215, 5608276.425 EPSG:32632 =

  837903.262, 5618484.501 EPSG:32631

why I get such a big difference: 5587608.614 instead of 5608276.425?
How does gdalwarp determine the pixel size and the extents?

EDIT:

The accepted answer shed light on my confusion. After I reprojected the lower right corner of the source raster instead of the lower left one, I got almost the same numbers:

x: 1093938.5119545625
y: 5618464.500556534
[668055.0297114132, 5587613.552738925]

x and y are the original coordinates in EPSG:32631


texts: GDAL creation option to set variable name to netCDF

texts:
Is there any way to set variable name in netCDF file with GDAL?
GDAL always sets the name Band1. I would like to set a different name. Still no success to find some -co (creation option) to overcome it.

E.g. command to rasterize a polygon

gdal_rasterize -a value -ot Int16 -te -180 -60 180 60 -tr 0.01 0.01 -of netCDF -co WRITE_LONLAT=YES -co TYPE_LONLAT=float -co FORMAT=NC4 -co COMPRESS=DEFLATE -l input input.shp result.nc

Possibly, is there any other simple way to do this?


texts: Lost ArcMap Add-in Toolbar Built with ArcObjects SDK in .NET and Problem with AssemblyCache folder and ArcGIS Cache Manager

texts:
Our organization manages an ArcMap 10.6.1 add-in written using the ArcObjects SDK for .NET Framework 4.7.2. We publish the add-in file to a server folder. We edited the registry to make ArcMap look for the add-in in the server folder.
REG ADD HKLM\Software\WOW6432Node\ESRI\Desktop10.6\Settings\AddInFolders /v "server_path_to_add_in_file" /t REG_SZ
The problem is that a handful, but not most, user machines lose the toolbar in ArcMap every time we publish an updated add-in to the server folder, and the user has to either reboot their machine or ctrl-alt-delete and find the ArcGIS Cache Manager and end that task to regain their toolbar.

I have read extensively about this problem:

https://community.esri.com/thread/121214#comment-625142
https://support.esri.com/en/technical-article/000012616
ArcMap Add-In Not Updating To Latest Version
https://community.esri.com/thread/117694

Interestingly, I see mixed results among users with their AssemblyCache folders (C:\Users\UserName\AppData\Local\ESRI\Desktop10.x\AssemblyCache). The add-in is associated with a GUID, and the folder name for the add-in in the AssemblyCache folder is that GUID. One user's AssemblyCache{GUID} folder updates to contain the latest dlls, and another user's folder is completely empty, always, before or after an update to the add-in.
I would like more information about this problem. Even better, I would like to know if there is an ESRI or Windows setting that can be changed to address this problem, or if that is not an option, is there a way to automate ending the ArcGIS Cache Manager task using ArcObjects inside the source code for our add-in?


texts: Comparing Tract-level data to State-level?

texts:
I'm using firearm trace data (available at the state-level) and attempting to compare it to socially disadvantaged areas (census tracts) by the way of regression analysis. Katz and Schnebly (2011) used population characteristics of tracts to measure the relationship between socially disadvantaged areas and gang members and/or gang crime.

I want to be able to measure whether the number of firearms being recovered in a state that are sourced to a different state are dependent on socially disadvantaged areas (certain neighborhoods in urban areas). Assuming that most firearms are being recovered in large urban areas due to higher levels of firearm violence in cities. Since my DV is only at the state-level, is there a meaningful way to generalize tract-level population characteristics to state-level?

I realize that one method could be summarizing these characteristics, but I wasn't sure if there would be any variation lost between multiple urban areas in states.

There is no way to disaggregate the state trace data to counties or tracts.

  Firearm Trace Data
  Recovered: (Factor) State where firearm was recovered and traced by law enforcement
  Year: (Integer) Year when trace request was initiated
  Source: (Factor) State where firearm was first purchased
  NumWeapons: (Numerical) Number of firearms that were recovered in state (Recovered) and successfully traced to a state (Source) during Year

The socially disadvantaged data consist of index scores (factor analysis using rates of poverty, single-families, 25 and over w/o HS degree, etc.) for each census tract in a state by year.


texts: How to specify the version of output SQLite database in ogr2ogr?

texts:
I'm having issues with a .sqlite database that I generate from PostGIS using ogr2ogr. The data appears to export into SQLite properly, and it loads in TileMill fine, but when I try to add any of the layers in QGIS it tells me:

Invalid layer. <layer name> is an invalid layer and cannot be loaded.

I was thinking that maybe my version of QGIS is running a slightly older version of SpatiaLite, so I decided to try exporting the data to a different version of SpatiaLite database. Is this an option in ogr2ogr?


texts: Recommended icon properties for QGIS plugin?

texts:
Starting with PyQGIS documentation - Developing Python Plugins:

icon - a file name or a relative path (relative to the base folder of the plugin’s compressed package) of a web friendly image (PNG, JPEG)

Sure icon is not the most important thing user cares about when needs to use your plugin to do some GIS processing etc. But it should look reasonably meaningful, and ideally nice, right? Any chance to use SVG? If not, and PNG/JPEG icon must be used as above - what should be the dimensions / color depth, any other recommendations?


texts:
I have a balloon callout ISymbol and I want to get a polygon object so I can see if the user clicks it.

I am trying to use the function QueryBoundary
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#//001w00000372000000

The function has four arguments:

public void QueryBoundary (
    int hDC,
    ITransformation displayTransform,
    IGeometry Geometry,  //what can I put in here?
    IPolygon boundary
);

One of my major hangups with the third argument is simply the fact that I am using this function to try to get the geometry... so what sense does it make to pass in the thing that I want?  Also since the function is called like this :

thissymbol.QueryBoundary(screenDisplay.hDC, displayTransformation, thisGeometry, polygon);

Presumably it doesn't need the geometry of "thissymbol" since "thissymbol" is the caller.


texts: QGIS 2.8.1 print composer won't export all OpenCycleMap base layer tiles

texts:
My set-up is Linux Mint 17.1, 32- bit system, QGIS 2.8.1-Wien.

I have a map with few layers including OpenCycleMap as a basemap.  In Print Composer the extent of the map is 32.5km x 22km square and the print page size is A1.

I am exporting from Print Composer to a large PNG image file of about 220Mb (container size) at a resolution of 200 DPI.

Everything exports to the image file fine, except OpenCycleMap. It does not appear to have enough time to download from the server.  There are therefore many missing background map tiles.


texts: Extract values from a field to a list in Python (arcpy)

texts:
I have a define function ( uniqueList ) in a module which sets the unique values from a list. I'm importing the function to a new script to set the unique values from a field 'COVER', after i've set the SearchCursor and when i run it, it doesn't work.

The define function :

def uniqueList(L):
    '''Find the set of unique entries in a list.'''
    S = set(L)
    # Copy S back to L
    U = list(S)
    return U

The script :

import sys, os
pathDir = os.path.abspath('C:\Gispy\Sample_scripts\Sample_scripts\ch16\exercises\listManage')
dirName = os.path.dirname(pathDir)
relativePath = 'listManage'
modulePath = os.path.join(dirName, relativePath)
sys.path.append(modulePath)
import listManager2
import arcpy
arcpy.env.workspace = 'C:/gispy/data/ch16/park.shp'
field = arcpy.da.SearchCursor('park.shp', ['COVER'])
value = listManager2.uniqueList(field)


texts:
Opengeo suite has provided a tutorial on how to manage LiDAR data using postgresql. I was trying to follow the tutorial in this link.

I am new to Linux and using Ubuntu 12.04. I have successfully installed Opengeo suite 4.1. I started off with the tutorial smoothly, installed point cloud extension and pdal and downloaded the example data. I only got stuck at the section "Loading LiDAR into the database". In this section, I followed the tutorial and implemented following codes successfully:

CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;

And also the following pdal commands run smoothly:

pdal info --input 20090429_42122c8225_ld_p23.laz --schema
pdal info --input 20090429_42122c8225_ld_p23.laz --metadata --xml

Now I faced the issue while implementing the following command:

pdal pipeline laz2pg.xml

From the terminal I got the error:

  mzaman@ubuntu:~/Downloads$ pdal pipeline laz2pg.xml
  Caught PDAL exception: unable to connect to database

This is to mention that, I have created a database in Postgresql called 'LiDAR' using pgAdmin. I kept another terminal open and run the following command using users 'postgres' and 'mzaman'.

mzaman@ubuntu:~/Downloads$ sudo -u mzaman psql LiDAR

psql (9.3.4)
Type "help" for help.
LiDAR=#

I also tried with the user 'postgres',

entermzaman@ubuntu:~/Downloads$ sudo -u postgres psql LiDAR

But the first terminal where I run pdal command, still shows the error:

  Caught PDAL exception: unable to connect to database.

My question is, how to bypass this error?


texts:
Is there a publicly available dataset somewhere that lists all country-adjacencies (pairs of countries that share a land border)?

Background: I want to programmatically access the list of countries that share a (land)border with country X.

Example: I programmatically want to access all of Luxembourgs neighboring countries, so having a dataset that contains these entries:

LU BE
LU DE
LU FR

Would allow my program to do useful stuff with the list [BE,DE,FR].


texts:
I found too many hits on the web to even see the right direction to go.

I have to provide a small proxy server serving requests for a line
of terrain elevation data between two coordinates.

It is a Linux system where GDAL is installed already. There are standard
raster elevation maps as inputs. The proxy needs to calculate about
1000 points/sec at runtime.

I'm not that much familiar with python but what I've read so far is that
GDAL/python seems the way to go. Is my feeling right to go this way?

thanks for feedback

Wolfgang R.


texts: Create a vector grid over a raster layer

texts:
I would like to create a vector grid over a raster layer, but with the grid being specifically only over the layer itself.

Right now, by using the vector grid option I'm getting a grid over a lot of unused space which is consuming considerable processing power and slowing down the work. Is it possible to create this grid but only over the parts of the raster that contain data?


texts: Python library for Geographically Weighted Regression (GWD)

texts:
The user @GrantD71 seems to have implemented a GWD in python through hard-coding using Numpy/Scipy. He asked for library/package to perform GWD but the best reply was given using an R package. I'm looking for a python library for GWD.

I looked in Scipy Docs and (almost) exhaustively googled it with no success. I found pygwr, but not sure if it's a completely developed and tested project. I even figured out how to install it and used it.

Does exists a wide-spread python library used for GWR? What is the most beneficial way (without the need of hard-coding) to perform it in python?


texts: Using sub-query in where statement of query in ArcGIS Map service?

texts:
I'm migrating my ArcGIS map services from version 10.0 to 10.2.   It uses ArcSDE to manage spatial data in a SQL Server database.  I was able to create the map service and everything seems OK. However if I add a sub-query in the WHERE statement of a query form in any map service, I receive "Unable to complete operation." error as the following image:

The same query was working in ArcGIS server 10.0. Also I can run the same query to select data, when I use the map document directly from ArcMap 10.2.
The ArcGIS Server 10.2, ArcSDE, and SQL Server are installed on a server which is internally accessible from the web server. The ArcGIS Web Adaptor is installed on the web server. Both servers have windows Server 2008 and IIS.


texts:
I'm looking for a way to find the max number of four different fields and write that in a new column. I have managed to use the field calculator to "+" the different fields into a new, but when I'm trying the max command I get lost. some of the fields are empty.

I have written this statement:

max(
(CASE WHEN  "mean" IS NOT NULL THEN "mean" ELSE 0 END),
(CASE WHEN  "_mean" IS NOT NULL THEN "_mean" ELSE 0 END),
(CASE WHEN  "mean_1" IS NOT NULL THEN "mean_1" ELSE 0 END),
(CASE WHEN  "max" IS NOT NULL THEN "max" ELSE 0 END),
)


texts: Storing XY coordinates in ArcGIS table as decimal numbers without losing zeros at the end

texts: I have created an ArcGIS table that contains coordinates as decimal numbers (float) with precision of 4 digits. The problem is that the zeros at the end of numbers in ArcGIS disappear, so I need a solution of how to keep/add zeros at the end of coordinates that have zeros at the end- for example, I need '64.12' to be stored like '64.1200'.


texts: How to download all addresses of some area from openstreet maps?

texts:
I need to know how to download all addresses of some area from openstreet maps?
After that I want to place all addresses in csv file.


texts: Select Duplicates in ArcMap Attribute Table

texts:
I am trying to select duplicates from a feature class in a file geodatabase using the ArcMap Attribute Table. The query I am using is:

Address_ID In (SELECT Address_ID FROM "FoodRoutes_20180615" GROUP BY Address_ID HAVING Count(*)>1 )

It is returning an invalid SQL statement error. Do the table or field names need to be in brackets? Anything I may be missing? I also tried the table name without " " and received the same error.


texts: How to add a layer on a RasterStack in a specific position?

texts:
I have the RasterStack called stack_raw. The dimensions are 517, 626, 323642, 46  (nrow, ncol, ncell, nlayers) .
I want to add an empty layer to a n position.
I´m trying this way, but it sounds a bit stupid:

#Creating the empty raster with the same dimensions:
a<-raster(nrow=517, ncol=626, ext=extent(stack_raw), crs=crs(stack_raw))
a[]<-1

#Separating the RasterStack in parts, in order to add the layer on that position:
stack_raw12<-subset(stack_raw, 1:2)
stack_raw346<-subset(stack_raw, 3:46)

#Adding the layer in the position of the middle:0
stack_raw12a<-addLayer(stack_raw12, a)

#Joining everybody:
stack_rawfull<-addLayer(stack_raw12a, stack_raw346)

Does anyone know a more intelligent way to do that?


texts: Connecting QGIS to remote Ubuntu PostGIS

texts:
I am trying to connect to a remote postgis 2.0 server (Ubuntu) from QGIS 2.0.1 (installed on winxp). I have been working with QGIS 1.7 / 1.8 on the same database for almost a year without any issue.

But I recently installed QGIS 2.0 and although everything seems to work fine, when I load a layer it has no rows at all, and when I try to exit QGIS it crashes. It raises a warning saying that:

  cursor 1 is lost

Could anybody give me some advice?

The log displays something like 20 PostGIS error messages, like this:

SQL: DECLARE qgisf0_11 BINARY CURSOR FOR SELECT st_asbinary(ST_Force2D("the_geom"),'NDR'),"id_column"FROM"schema"."table"WHERE"the_geom"&&st_makeenvelope(41.5457,35.544,36.784254,44.15748)
Resultado:7(ERROR: function st_force2d(geometry) does not exist
LINE1: ...RE qgis0_11 BINARY CURSOR FOR SELECT st_asbinary(ST_Force2D...
HINT: No function matches the given name and argument types. You might need to add explicit type casts.

The actual name of the function in PostGIS 2.0, and also in postgis 1.5 is st_force_2d(geometry)


texts: World_Eckert_IV to Latitude/Longitude projection is incorrect

texts:
I'm using pyproj to transform between coordinate systems.
Here is my .prj file:

PROJCS["World_Eckert_IV",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Eckert_IV"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]]

where I call

project = partial(pyproj.transform,
                  pyproj.Proj(proj="eck4"),
                  pyproj.Proj(init='epsg:4326'))

output["Latitude"] =\
    dataframe['geometry'].apply(lambda x: transform(project, x)).y
output["Longitude"] =\
    dataframe['geometry'].apply(lambda x: transform(project, x)).x

However, the latitude/longitude coordinates that I get back do not point to where they should. Does anybody have any ideas?


texts: How to specify (cross-browser) mouse pointer sets on leaflet control?

texts:
Leaflet custom controls might be interactive, when the user can use the mouse on them to perform an action, or non-interactive, when they are used only to show information. This should affect the mouse pointer used for elements of the control (and possibly sub-elements of it for partly interactive controls).

What is the correct way to set those pointers for the control container and its sub-elements, if any (e.g. well defined css classes or else)?


texts: which map projection should be used for global map with geography and human demographics

texts: I am working on a project to try to identify the most threatened populations (on a global scale) of having a limited access to freshwater. I want to compare those most at risk currently and then add in climate change projections to see how it is expected to change. I am trying to justify what projection system to use. With such a global scale what is the best option, if focusing on freshwater geography as well as human demographics? I am a little overwhelmed, but looking at mostly WGS 1984 coordinate systems. How would you justify a certain projection for such a broad analysis area and concentrations.


texts: Can't Find Option to Create Feature Layer

texts:
I have created an ArcGIS online account but I can't find option to create feature layer.

I tried to create developers account, but yet didn't find a way to create feature layer service.

I have added map from CSV file, but don't find publish button so that it can be hosted using feature layer.

How do I host my service as feature layer?


texts: Calculating Area of Intersecting Polygons and Buffers

texts: SO, in QGIS I have several points with 1km buffers around them; some of the buffers overlap. Within the buffers are a number of ponds, which I have a shapefile for. I want to calculate the area of pond within the buffers, but if one pond is covered by overlapping buffers, I do not want the area counted multiple times, which is the issue I'm dealing with right now. Can this be done? Creating an intersect shapefile then adding geometry columns counts overlapped ponds multiple times.


texts: Calculating Accuracy Assessment in QGIS?

texts:
Is there a way to calculate an accuracy assessment in QGIS?

I have around 300 points. At each point I extracted the value of the pixel. I compared the extracted value with the value which should be there. Now I have a table with two columns: One is the extracted value from the raster, one is the value which I typed in manually (the correct value).

How would I calculate the users accuracy and Kappa with that? Is there a tool to do that?


texts: Exporting QGIS Field Calculator user defined function

texts:
I have created many Python functions in QGIS Field Calculator Function Editor. How can I export the user defined functions to use in other computers or share to other people?
My QGIS is 3.12.1 on MacOS 10.15.7


texts: Modifying null values in ArcPy for ArcGIS Pro?

texts:
How can I say in ArcGIS Pro if the field (String) is Null, then replace it with "No"?

I was trying this:

database_field.replace("Null","No")

No error occurred, but the null values did not change...


texts: Recommendations about graphic interface for ArcPy/Python script

texts:
I have done some scripting using ArcPy and Python, and I want to make some interface to select some shapefile. What are the best graphic tools for Python to make this?

I don't need high detailed graphics.


texts: A missing feature class in the list when adding classes to a new network dataset

texts:
I´m having a problem here when creating a new Network Dataset. Firstly I select the feature classes from the list that will participate in the network dataset. I have a feature dataset under a geodatabase file, which stores different road elements (major roads, private roads, biking and walking ways etc.) and lastly a feature class, where the information of all previously listed different road layers are in. That´s why this feature class layer is  important. I have made some modifications to this layer: I added fields in the attribute table that for example show the travel time by bicycling.

By some reason the New Network Dataset -wizard doesn´t find this feature class, that I´ve modified. Can you tell me why? Might it be the reason that I´ve changed it a bit and somehow that has also modified the way the feature dataset is being built? Do I have build it again somehow?


texts: UpdateLayer function in arcpy - runs, does not update layer

texts:
I'm very new to ArcPy and programming in general. I'm trying to use UpdateLayer  to add hatching to a selected layer from a template layer I've created.

The script runs (no errors), but it doesn't update my selected layer:

import arcpy

mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, 'Stationing')[0]
sourceLayer = arcpy.mapping.Layer(r'H:\Python\Data\Template.lyr')
updateLayer = arcpy.mapping.ListLayers(mxd, 'IAmRoute.lyr', df)[0]
arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, False)

mxd.saveACopy(r'H:\Python\Stationing2.mxd')
del mxd, sourceLayer

Am using ArcGIS Desktop 10.5.1.

When I save a copy of the mxd, I end up with two layers named Template. One appears to be the original Template layer, the second appears to be my original IAmRoute layer with its original properties.

Running from Arcmap w/mxd Stationing open.

edit: changed MapDocument to 'CURRENT'


texts: Dijkstra/Finding shortest path gives back a tree with only -1 in it

texts:
So I'm trying to implement a routing algorithm in my QGIS 3 plugin that I'm writing using Python.

I think I could resolve all the problems concerning the break between QGIS 2 and 3, but my resulting tree gives me for all nodes in it just a -1 value which is for not connected nodes.

If I have for example a project that looks like this:

Where the red Square is my start point and the yellow circle my end point I get back a tree that looks like this [-1,-1,-1]

the code I use is the following:

First I select the Line layer and the startPt and endPt

director = QgsVectorLayerDirector(layer, -1, '', '', '', 3)

strategy = QgsNetworkDistanceStrategy()

director.addStrategy(strategy)

builder = QgsGraphBuilder(crs) #crs = project.crs() which is in this case WGS84

tiedPoints = director.makeGraph(builder, [startPt, endPt])
graph = builder.graph()

startVertex = graph.findVertex(startPt)

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startVertex,0)

I just don't know where it goes wrong.


texts: Relative path not working in QGIS with Google Drive?

texts:
I have a qgis file, referring to a few shp files, saved on a google drive folder which is synced across two computers using different OS, one is Mac OS X, the other Windows 8. So, the absolute path of this qgis file on different operating systems is different. But the location of this qgis project file in relation to the shp files it refers to is the same.

E.g:
On the Mac OS X computer, the locations of the project file and the shp file are:

/Users/David/GoogleDrive/GIS/Project/project_01.qgs

and

/Users/David/GoogleDrive/GIS/Project/data/road.shp

On the Windows 8 computer, the locations of the project file and the shp file are:

C:\Users\David\GoogleDrive\GIS\Project\project_01.qgs

and

C:\Users\David\GoogleDrive\GIS\Project\data\road.shp

I setup this qgis project file as using relative path, following the suggestion from the following post:
relative path

This project works file in Mac OS X, even if I move the qgis file and the associated shp files to different directory.

However, I got a window popping out when opening it Windows 8, saying "Handle bad layers" and asking me to re-specify the locations of all the shp files again.

Isn't the relative path supposed to deal with this kind of situation?

May I ask how to solve this problem related to qgis project file saved in Google Drive folder that is synced across different operation systems?


texts:
Is there any way to have multiple text formats in a single qgis label?  I am making a map that shows mountain summits.  Many of these summits have more than one name and I would like to show each name in a label and have each name rendered with different formatting to show what type the name is.

For example, summits, in the US, often have an "official" name from the GNIS database, a name as shown on a 1:24000 USGS topo map (often just the elevation of the summit), and a local nickname.  I would like to be able to have a label that renders each name in a different type font and/or font size.


texts:
I downloaded from the USGS Earth Explorer a MODIS Landsurface Reflectance V6 MODIS MOD09Q1 V6 image to calculate an NDVI. I'm trying to open it on ArcGIS Pro, however when I import the file, the Subdataset Selection dialog box does not open and the hdf files appears as a time serieS.
Does anyone know how to properly open this file or how can I open the Subdataset Selection dialog box?


texts: Calculate nearest distance in ArcView, XY coordinates

texts:
I am really lost and struggling because I have never used an Arc system before.

I have two datasets. One with house locations, and the other with investment locations- each locating by XY coordinate.

The room where I conduct my analysis only has ArcView, so I must use that.

I want to calculate the distance between a house and the nearest investment (calculate the distance between two XY coordinates).

Is this possible with ArcView?

I also want to join the variables in the investment dataset to the housing dataset, after I have calculated the nearest distance. Is that possible?

I am sorry if I am not explaining myself correctly- I am just very confused, if it's not stata, then I don't know how to use it!


texts: Opening .mxd files in Python not via ArcPy?

texts:
I am an old user of Python, but relatively new to GIS. I have a .mxd file and would like to read it into Python.

Is arcpy the only way to do so?

I tried with pyshp but it does not accept this file extension.

Is there perhaps another package for python3 for working with data in this format?


texts: Integrate variable in GDAL raster calculator in script editor QGIS 2.8

texts:
I'm writing a script that will iterate rasters and do some calculations on each raster. I would like the formula in the GDAL raster calculator to include a variable, which I read from a txt file. I don't want the user to insert the variable but to integrate it in the script.

The code of the raster calculator I use is this:

outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', rlayer,'1',None,'1',None,'1',None,'1',None,'1',None,'1','10000*(A*0.00002-0.1)/sin(65)*(A>0.001)','0',4,None,outputName)

I need the number inside of the sin, in this example 65, to be a variable, so every time I'll have a different number.

Is that possible to integrate a variable in this syntax?


texts: Exporting a grid to KMZ with legend (QGIS 3)

texts:
I am stuck trying to export a grid with values to a KMZ file with a legend that can be viewed in Google Earth Pro.

This is the grid I'm trying to export (it is in EPSG 32756)

This is how I need the KMZ to be displayed and structured in GEP:

And this is how the attribute table is constructed:

I've tried exporting to KML using polygon grometry but even if I convert the CRS to EPSG 4326 for GEP first, it ends up on the other side of the world, or on top of the north pole, and with pretty much no data in it. Can anyone give me some ideas of how I should be exporting this data?


texts: Calculating distance from point to feature returned from layer query using ArcGIS API for JavaScript

texts:
I have a project where a user clicks on a map and I filter to show all features from a layer that are within a certain distance of the point clicked.  That part of things works fine.  What is not working are my attempts to calculate the distance from the point clicked to the queried feature location.  I've isolated the issue to my geometry service:

export class GeometryService {
public geodesicUtils;
public Point;

constructor() {
    this.initialize
}

private async initialize() {
    try{
        [
            this.geodesicUtils,
            this.Point
        ] = await loadModules([
            GeodesicUtilsModule,
            PointModule
        ]);
    } catch(err) {
        console.log(err);
    }
}

public getDistance(from: {x: number, y: number}, to: {x: number, y: number}, unit?: string): esri.GeodesicDistanceResult {
    const fromPoint = new this.Point(from);
    const toPoint = new this.Point(to);
    const result: esri.GeodesicDistanceResult = this.geodesicUtils.geodesicDistance(fromPoint, toPoint, unit);
    return result;
}

Essentially the app just stops when instantiating a new point.  When I break point through it just goes off into some obfuscated dojo nightmare code and never returns a result.  No exception is thrown, no error displayed it just disappears into purgatory never to be heard from again. Below is the code that initiates the call to the service:

this.queryResults = results.features.map((feature: any) => {
    const deviceFeature: IElectricDeviceFeature = this.electricDeviceFeatures.find((electricFeature: IElectricDeviceFeature) => electricFeature.value == feature.attributes.assetgroup);
    const toPoint: {x: number, y: number} = {
      x: feature.geometry.longitude,
      y: feature.geometry.latitude,
    }
    const calculatedDistance: esri.GeodesicDistanceResult = this.geometryService.getDistance(fromPoint, toPoint, this.searchState.units);
    return {
      attributes: feature.attributes,
      geometry: feature.geometry,
      symbol: deviceFeature.symbol,
      label: deviceFeature.label,
      distance: calculatedDistance.distance,
      distanceUnits: this.searchState.units,
    };
  });

This is an angular 8 app using the arc gis 4.14 javascript api.  It also uses some ngrx for statemanagement although that is not relevant here.


texts: MapServer Mapfile and OpenLayers: using SIZEUNITS. Labels shrinking on 1x zoom and points disappearing upon zooming in

texts:
I am trying to set the symbol radius to 100 miles. This seems ok (note, not including the whole mapfile, just what's at hand here):

LAYER
  NAME POIs
   TOLERANCE 50
   TOLERANCEUNITS miles
   LABELITEM "name"
   SIZEUNITS miles
  CLASS
    STYLE
       SYMBOL 'Airport'
       COLOR 61 245 0
       SIZE 100
       MINWIDTH 100
       MAXWIDTH 100
    END
    LABEL
      SIZE 9
    END
  END

I'd like to keep the label size at 9. If I am NOT using

SIZEUNITS miles

the label at size 9 looks fine. But using SIZEUNITS miles makes it much much smaller.

After reading here: http://lists.osgeo.org/pipermail/mapserver-users/2012-February/071582.html it seems you can set a maxwidth and minwidth within your STYLE to accompany for having the symbol width in pixels but the radius in miles (correct me if I am wrong).

Also:

When I zoom in on a handful of points, certain points and labels start to disappear. I know I could set FORCE ON and PARTIALS TRUE but this doesn't look good.

So what I would like to understand:

1) why are the labels becoming much smaller with SIZEUNITS miles on?

and

2) why are certain points disappearing when zooming in (and some only showing up when zooming in, not initially present) when there seems to be ample room on the map for them to display without FORCE ON? . Right click image link into new tab to see better view.


texts: Why is short-wave infrared portion of electromagnetic spectrum so sensitive to fire?

texts:
After attending a remote sensing workshop and using Copernicus data to monitor natural disaster scenarios, my question could not be fully answered by the instructor. All that was said, was that is in no way associated to thermal IR.

What property does fire have, that makes is visible with shortwave infrared?


texts: Splitting Polygon using Polygon in ArcGIS for Desktop?

texts:
I want to split the polygon layer Settlement using another polygon layer District boundary.

How can I do this in ArcGIS ?

There are 64 districts adjacent to each other. Different part of a Settlement might be in 2/3 adjacent districts, which I need to split.


texts: Capturing multiple screenshots of image at extents from point features in ArcGIS Desktop?

texts:
I've got a map with ~2000 villages on it & have a layer of points that corresponds to them. I'd now like to enter some information about them (the names). It would be a huge pain to do this manually by scrolling around the layer in Arcmap.

One idea would be to programmatically take screenshots in a little box around each point. I could then do data entry based on the screenshots pretty easily.

Is it possible to do this with ArcGIS Desktop?

I have a .tif image I've georeferenced and placed points on, sample as image below. I'd now like to do some data entry for each of the points. The simplest thing would be to take a screenshot of a certain area around each point and scroll through them.

The ideal final product would be:

2000 individual image files
Each a screenshot centered around a point in my file with a certain zoom factor applied
The below example would generate 8 image files, each centered around a red dot with a certain zoom (example output below)


texts: how to remove previous marker automatically before adding new marker

texts:
I am using ui-leaflet on my ionic 1 mobile app which is the same as angular-leaflet-directive shown in the plunker link.

I am getting the user location and adding it to the marker and circle in the map, the issue i am having is that has the user moves their location is updated and therefore adding a new marking which then spans the maps which is not what I want.
What I am trying to achieve is that has a new marker is added or before a new marker is added, It should remove the old marker. How can I do this please? Below is my code and what I have tried

var app = angular.module("myApp", ['leaflet-directive'])

.controller('myCtrl', myCtrl);
myCtrl.$inject = ['$scope', 'leafletData'];

function myCtrl($scope, leafletData) {
  var location = {
    lat: -33.8830,
    lng: 151.2166
  };

  var mainMarker = {
    lat: location.lat,
    lng: location.lng,
    focus: true,
    draggable: false
  };

  var vm = angular.extend(this, {
    center: {
      lat: location.lat,
      lng: location.lng,
      zoom: 17
    },
    markers: {
      mainMarker: angular.copy(mainMarker)
    },
    defaults: {
      zoomControl: false
    },
    tiles: {
      url: 'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png'
    }
  });

  function onLocationFound(e) {
    var radius = e.accuracy / 2;
    var marker;

    leafletData.getMap().then(function(map) {
      /* does not work in removing the old marker */

      if (marker) { // check
        map.removeControl(marker); // remove
      }

      var options = {
        icon: {
          icon: 'fa-certificate',
          prefix: 'fa',
          shape: 'square',
          markerColor: 'green'
        }
      };

      marker = L.marker(e.latlng).addTo(map);
      //.bindPopup("You are within " + radius + " meters from this point").openPopup();

      L.circle(e.latlng, radius).addTo(map);

    });
  }

  function onLocationError(e) {
    alert(e.message, e.code);
  }

  leafletData.getMap().then(function(map) {
    L.easyButton({
      position: 'bottomright',
      states: [{
        icon: 'fa fa-location-arrow fa-lg',
        onClick: function(btn, map) {
          map.locate({
            setView: true,
            maxZoom: 16,
            watch: true,
            enableHighAccuracy: true
          });
          map.on('locationfound', onLocationFound); //see if it is possible to pass the map has a perimetrr
          map.on('locationerror ', onLocationError);
        }
      }]
    }).addTo(map);
  });
}


texts: Sharing QGIS projects with Geopackage connections

texts: I've built a QGIS template project that I'm using for repetitive municipal projects that I'm sharing with other users on the local network, the layers are organized and stylized within geopackage files.  I've just learned that when other users open the QGIS project, they don't see the layers because they're not connected to the geopackage files - not a hard problem to solve since they can just make the connection and move on, but I feel like there must be a way to have the qgis project automatically make that connection so it doesn't have to be done manually for each user?


texts:
I want to delete about 25 datasets that are empty from my gdb file which consists 30 datasets. I am invoking pPthon window and trying to execute this script:
import arcpy
from arcpy import env
env.workspace =r"D:\test\E44M11_HYD_UTM44.gdb"

data_sets = arcpy.ListDatasets()

for item in data_sets:
    fcList = arcpy.ListFeatureClasses('', '', item)
    for item in fcList:
        fcLength = arcpy.GetCount_management(item)
        if int(fcLength.getOutput(0)) == 0:
            arcpy.Delete_management(item)

But just after 2 layers deletion it prompts this error message.

How can I fix this issue?