texts: Trigger on PostGIS line layer occurs twice

texts:
Using this tutorial, I've established a SQL trigger in order to track changes on a PostGIS layer. It works great for a point layer. However, using the script below, it causes problems with line layer, as all of the triggers occur twice.

I first realized this because the history table inserts lines for the 'INSERT' and 'UPDATE' logic twice (not for 'DELETE', surprisingly). The notices print to the console twice for 'INSERT' 'UPDATE' and 'DELETE'.

For example:

When creating a line in the base table, the history table creates two objects with the same attributes except for the hid (history table id).
When updating the attribute of an existing line, in the history table, three objects are created:
1) One with etat = 'MODIFICATION COURANTE' (resembles the existing feature in the base table, including the updated attribute)
2) Two with etat = 'MODIFICATION ARCHIVEE' -- one with the old attribute value and one with the new attribute value. At this point, I am expecting only one object containing the old attribute. The object containing the new attribute is out if place.

--CREATE TRIGGER
CREATE or REPLACE FUNCTION test.test_track_history_tracker() RETURNS trigger AS
$new_test_track_history_tracker$
            BEGIN
            -- INSERT
            IF (TG_OP = 'INSERT') THEN
                           INSERT INTO test."Conduite_test_history"
                                           (diametre, type, materiau, origininfo, dtreno, dtpose, reparateur, dt_dmd, id_sig, gid_org, created, created_by, modified, etat, geom)
                           VALUES
                                           (NEW.diametre, NEW.type, NEW.materiau, NEW.origininfo, NEW.dtreno, NEW.dtpose, NEW.reparateur, NEW.dt_dmd, NEW.id_sig, NEW.gid, current_timestamp,
                                           current_user, FALSE, 'CREATION', NEW.geom);
                           raise notice 'Insert happened';
                           RETURN NEW;

            -- UPDATE
            ELSEIF (TG_OP = 'UPDATE') THEN
                           UPDATE test."Conduite_test_history"
                                           SET deleted = current_timestamp, deleted_by = current_user, modified = TRUE, etat = 'MODIFICATION ARCHIVEE'
                                           WHERE deleted IS NULL and gid_org = OLD.gid;
                           INSERT INTO test."Conduite_test_history"
                                           (diametre, type, materiau, origininfo, dtreno, dtpose, reparateur, dt_dmd, id_sig, gid_org, created, created_by, modified, etat, geom)
                           VALUES
                                           (NEW.diametre, NEW.type, NEW.materiau, NEW.origininfo, NEW.dtreno, NEW.dtpose, NEW.reparateur, NEW.dt_dmd, NEW.id_sig, NEW.gid, current_timestamp,
                                           current_user, FALSE, 'MODIFICATION COURANTE', NEW.geom);
                           raise notice 'Update happened';
                           RETURN NEW;

            -- DELETE
            ELSEIF (TG_OP = 'DELETE') THEN
                           UPDATE test."Conduite_test_history"
                                           SET deleted = current_timestamp, deleted_by = current_user, etat = 'SUPPRESSION'
                                           WHERE deleted is NULL and gid_org = OLD.gid;
                           raise notice 'Delete happened';
                           RETURN NULL;
            END IF;
END;

$new_test_track_history_tracker$ LANGUAGE plpgsql;

--ADD TRIGGER
DROP TRIGGER IF EXISTS trg_test_track_history_tracker ON test."Conduite_test";
CREATE TRIGGER trg_test_track_history_tracker AFTER INSERT OR UPDATE OR DELETE ON
test."Conduite_test"
            FOR EACH ROW EXECUTE PROCEDURE test.test_track_history_tracker();


texts: ArcObjects: How to get the midpoint of a set of polylines?

texts: How can I find the midpoint (along the curve) of a set of polylines using ArcObjects?


texts: QGIS 3.2 - Poor conversion from high-resolution vector to raster

texts:
I have a shapefile with polygons with several sizes, and I want to convert it to a raster with 30mx30m pixel resolution.

Although I was able to obtain the raster with the pretended resolution (using QGIS 3.2,Raster > Conversion > Rasterize), the file I obtained ignores a lot of several small polygons; that is, the polygons are not converted in one or more pixels.

Additionally, there are small polygons converted in a pixel but other larger polygons are not converted. See this example below:

In this case, the pixel is 30x30m, A is ~111 and B is ~1642 square meters.

Questions:

1) Does anyone know how this conversion works, i.e., what are the criteria for converting a polygon in a raster?
Using the example in the picture, why was A converted in a pixel and B was not?

2) Is there a way to calibrate/choose the size of the polygons to be considered for conversion?


texts:
I'm trying to add a GUID field programmatically to an existing FeatureClass using ArcObjects.

Here is my code:

ISchemaLock schemaLock = (ISchemaLock)ipFC;

try
{
    // A try block is necessary, as an exclusive lock may not be available.
    schemaLock.ChangeSchemaLock(esriSchemaLock.esriExclusiveSchemaLock);

    //Ajout d'un champ GUID
    IFieldsEdit pNewFieldsEdit = (IFieldsEdit)ipFC.Fields;
    IField guidField = new FieldClass();
    IFieldEdit guidFieldEdit = (IFieldEdit)guidField;
    guidFieldEdit.Name_2 = "OldGlobalID";
    guidFieldEdit.Type_2 = esriFieldType.esriFieldTypeGUID;
    guidFieldEdit.DefaultValue_2 = null;
    guidFieldEdit.IsNullable_2 = true;

    pNewFieldsEdit.AddField(guidField);
}
catch (Exception exc)
{
    // Handle this in a way appropriate to your application.
Console.WriteLine(exc.Message);
}
finally
{
    // Set the lock to shared, whether or not an error occurred.
    schemaLock.ChangeSchemaLock(esriSchemaLock.esriSharedSchemaLock);
}

Everything seems to work (no error) and it's based on an documentation example. Unfortunately, when I click on the feature class in ArcCatalog and open the table preview, I have this message showing:

Also, the table open in the preview, but it's all empty (header and row) with 308 records(the correct number of row in the feature class). What am I doing wrong?

EDIT

If I check the properties of my feature class, the new field seems to be there:

But the preview looks like this right after the error message...:

Any ideas?

EDIT 2

I discovered that when I restart ArcCatalog and go back in the properties, the new field is gone and there is no problem opening the table preview, but the new field is nowhere to be found. Maybe there is a function to "Save" everything after the field has been added?


texts: Using Engineering symbols (*.svg) in QGIS?

texts:
QGIS looks straight forward except for the symbols or icons.

I have created SVG symbols and when I place my pointer to the obvious blank areas in the symbol library a pop-up says that my symbol is there, yet nothing shows in the icon library.

I insert the "ghost" svg into the map and Nothing shows on the map.

Can anyone tell me how to get an SVG file to work in QGIS?


texts:
I am working with a RasterStack, from which I only need the data for a certain region, defined in a SpatialPolygonsDataFrame. I used the following code which works well:
test.ex <- crop(dat.stack, extent(myregion))
test <- mask(test.ex, myregion)

Is it however possible to add a "buffer" around my area of interest, so that not only the grid cells of the RasterStack covered by the polgyon of my region but also one additional grid cell are added around the outline of the resulting raster?


texts: Connect Points using line different table in QGIS3

texts:
I have two tables.

Table1 is having many point objects (Sample given below)
Table2 shows the connectivity between the point object(using a link) and link details

My requirement is to have a new Line layer connecting point objects from Table1 as per the Link in Table2.

example : I need to have one line connecting A1 to A3 with a name of "A1_A3", another Line connecting A5 to A4 name of the line to be A5_A4)

I tried with few plugins like connectpoints & points2one... But no luck.


texts: d3.geo.path() to draw a path from GIS coordinates

texts:
Trying to take GIS data of the zip codes of the Bronx and turn it into a map using d3.js. But I am having trouble rendering. Could you show me proper use of d3.geo.path() ?

[ [ -73.84463693999993, 40.90475939500004 ], [ -73.84443733899991, 40.90423928800004 ], [ -73.84443555299993, 40.904234630000076 ], [ -73.84332154899994, 40.90397722300008 ],...]

These numbesr are taken from NYC zip codes map: https://nycopendata.socrata.com/Social-Services/Zip-Codes-Map/zsjh-u7ve

In fact, this question goes for any list of GIS coordinates.  How do I turn it into a path?  I could not convert the documentation this into code. https://github.com/mbostock/d3/wiki/Geo-Paths


texts: How to display X-Y data with different colors using arcmap

texts: I have a list of lat/lon coordinates for data points that I want to put on a map.  I want to have the size and/or color of the points to change based on additional columns on the excel sheet with the coordinates (kind of like a bubble chart).  Is it possible to do this in GIS? I am using arcmap 10.2.


texts: Using ECQL filter on user form using GeoServer and OpenLayers

texts:
I know this is more of a JS question (skills of which I am inadequately inept) but it is using OpenLayers and GeoServer.

I am trying to apply a filter with a user form on a Layer, but am not sure why it is not working.  If I manually write the code into the CQL_FILTER, it works, but it won't work with the UpdateFilter button.

function init_map() {

    var layerExtent = new OpenLayers.Bounds( -9816656, 5296800, -9772016, 5345759);

    map = new OpenLayers.Map('map',
            {'restrictedExtent': layerExtent});
    basemap = new OpenLayers.Layer.OSM("Simple OSM Map");
    map.addLayer(basemap);
    markers = new OpenLayers.Layer.Markers("Markers");
    map.addLayer(markers);
    map.setCenter(
        new OpenLayers.LonLat(-71.147, 42.472).transform(
            new OpenLayers.Projection("EPSG:4326"),
            map.getProjectionObject()
        ), 12
    );

    var tiled = new OpenLayers.Layer.WMS(
                "Verizon Cell Towers", "http://mpdflash:7070/geoserver/Milwaukee/wms",
                {
                    LAYERS: 'Milwaukee:Verizon Cell Towers',
                    tiled: true,
                    transparent: true
                },
                {
                    displayOutsideMaxExtent: true,
                    isBaseLayer: false,
                    yx : {'EPSG:6948' : false}
                }
            );

    map.addLayer(tiled);

    tiled.mergeNewParams({'CQL_FILTER': "filter"})
}

function updateFilter(){
        var filter = document.getElementById('filter').value;
    }

function resetFilter() {
        document.getElementById('filter').value = "";
        updateFilter();
    }

and

            <input type="text" size="80" id="filter"/>


texts: Unable to edit file geodatabase: How do I fix it?

texts:
I have different geodatabase for different project. When I start projectA.mxd which connects to projectA.gdb, I cannot create and edit feature class inside. However, when I start projectB.mxd, I can create and edit projectA.gdb.

Why???? How can I fix it?

I am using ArcMap 10.1.


texts:
I want to iterate a buffer from a list of distances I have, but I can´t find the way to do it. I´m a beginner on Python.

import processing
x=[100, 200, 300, 400, 500]
processing.run("native:buffer",\
{'INPUT':'C:/Path/points.shp',\
'DISTANCE':x,\
'SEGMENTS':5,\
'END_CAP_STYLE':0,\
'JOIN_STYLE':0,\
'MITER_LIMIT':2,\
'DISSOLVE':False,\
'OUTPUT':'"C:/Path/Buffer" + str(x) +".shp"'})


texts: I have created a polygon from grass gis and found small triangle like structure inside the polygon (sometimes at the edge to outside), that I want to remove/fill. I tried with geometry check to fix but not succeed. Is there any efficient way to fix that.


texts: Efficiently checking for feature overlaps with ArcPy?

texts:
In an effort to reduce processing time, I am looking for ways to pre-process data within my script. I have ~18,000 featureclasses with ~500 features each (all simple circle polygons) from which I want to eliminate coastal land area (another simplified polygon) using the Erase function.
My idea for reducing overall processing time is to check for overlap between each feature class and the coastline polygon. If there is overlap, the script should erase the land from the features. If there is no overlap, the script should move on, and look at the next feature class.
Specific questions:

What is the most efficient, Pythonic way to say "if there is overlap, do XX?" I cannot find documentation to do something similar to "if Overlap == True" besides the method:
MakeFeatureLayer-->SelectFeatureByLocation-->CopyFeature-->GetRowCount--> if RowCount>0

Would this actually reduce processing time or am I crazy?

I'm working locally on my machine, running the script in IDLE, with ArcGIS 10.1 (with all licenses) on a 32-bit Windows 7, 3.33GHz Intel DuoCore, 4GB RAM computer.


texts: How do I install and run SEXTANTE for QGIS 2.0?

texts: I just upgraded to QGIS 2.0, and I have been unable to locate/use SEXTANTE tools.  A search with the plugin installer does not turn up anything for "sextante."  Any idea what is going on?


texts: Finding distance from river outlet to points, following river, using ArcGIS Desktop?

texts:
I have a line (the river (in blue)) and a shapefile of points (plants colonies (in orange)) and I need to find the distance between the river outlet and many points along the river instead of "as the crow flies". I want to determine the order of the colonies, downstream to upstream, and their exact position in meters along the river. I'm using ArcMap 10.4 and I have all extensions needed (Network analyst, Spatial analyst, etc.).

I read many forum topics about that and I think locating features along routes is what I need (Do you agree ?). I created a "Network Dataset" from my river line to make a "new route" in the Network Analyst toolbar, which is supposed to make me a "route analysis layer" that I could use to locate my colonies along the river. But the attribute table of the route analysis layer is empty and so the geoprocess gives nothing...

Do you know what I'm doing wrong or what should I do to calculate the distances I need ?


texts: Changing field formats in ArcGIS geodatabase?

texts:
I am working with data in an ArcGIS 10.1 geodatabase.  I wanted to complete the simple task of taking a ten-digit ID field that is currently in String/Text format and create an identical field in Long/Numeric format.

However, Esri's resource pages on tables (http://resources.arcgis.com/en/help/main/10.1/index.html#/Common_tables_and_attributes_tasks/005s00000005000000/) do not discuss how to achieve this.  Also, the following GIS/SE thread does not help me either: Is there a way to convert multiple columns at once from string to text, or to other data formats? and the post mentioned in the comment there refers only to QGIS.

In the meantime, I created a new table with the extra field, joined it to the original shapefile and will reimport this shapefile into the gdb, adjusting the new field format as needed.

Is there an easier way to complete this task?


texts: Calculating gain and loss area of vegetation land from two lines representing different years using ArcGIS for Desktop?

texts:
How can I calculate a gain and loss area of vegetation land having two different lines that represents vegetation in different years.

Edit:

There are two polylines, one represents the edge of vegetation of a forest in 2010. The other one the edge in 2007. I need to create a polygon of areas where there is a gain or loss in land area.


texts: Changing values of Leaflet MiniMap after creation

texts:
I want to change the values of Leaflet MiniMap object after it is already created.

e.g. miniMap.options.width:

var miniMap = new L.Control.MiniMap(
  tileLayer01, { centerFixed: new L.LatLng(50.0, 20.0) }
  ).addTo(map01);

console.log("miniMap.options.width: " + miniMap.options.width);
// prints 150 (the default value the minimap was initialized with)

miniMap.options.width = 80;

console.log("miniMap.options.width: " + miniMap.options.width);
// prints 80 (the newly set value)

However, that has no effect, the minimap is not rendered accordingly but stays unchanged.

How can I change MiniMap values after creation so that it affects the minimap?


texts: Is there a standard measure of Remoteness?

texts: I have a number of sites around the globe, and they are not uniformly distributed. I would like to use some measure of remoteness, in order to compare sites that are more or less remote. I guess something like 1/sum(distances to other sites) would work, but I would first like to know whether there are any other standard metrics for measuring the remoteness of individual sites against a large population?


texts: Joining attributes of polygons and points by location in QGIS

texts:
I have two shapefiles....

points containing a number field 'response_time'
a polygon grid of square polygons containing fields including 'PNTCNT' with the number of points per polygon, calculated already using the tool 'points in polygon'

There are between 0 and 600 points within the bounds of each polygon in the grid.

I want to add fields to the polygon grid for the mean, min and max response times within each polygon, in addition to the PNTCNT field already there.

I have tried the 'Join attributes by location' tool, making the polygon layer the target and the points the join layer, and vice versa. I've selected 'take summary of interesting features'.

But I can't get the result I want. All I end up with is the polygon grid again but with an extra field 'COUNT' which replicates (or duplicates) the 'PNTCNT' field. No min, max and mean of values from the points layer.

What am I doing wrong? Is there a better way to do this? I have tried this using QGIS 1.7 and 1.8.


texts: Using PostGIS functions in ArcGIS Desktop for Layer Definition Query?

texts:
I am connecting ArcGIS Desktop 10.2 to Postgis 2.0.3 with PostgreSQL 9.2.4. No SDE. I can load data with no problems and access the data in ArcMap. However I can't seem to use PostGIS functions within a layer definition query. For example:

ST_DWithin(location, ST_GeomFromText('POINT(-97.5 35.5)', 4326), 0.01);

Which works perfectly from PostgreSQL returns the error:

  One or more layers failed to draw:
  Report:  Failed to parse the where clause.

This is strange, as if the exact same query is used to define the underlying source query for the layer it works fine. Am I missing something here. I used to have this working for Oracle no problem. What doesn't it like about PostgreSQL's query?


texts: Creating new Half-Inch Bartholomew-Style Maps in QGIS 2.18?

texts:
I have a full set of the old half-Inch Scale John Bartholomew maps produced up to the early 1970s (1972 - 1974) of Great Britain, and I'd like to produce new mapping in the same style, but with modern Sans Serif settlement labels and fully up-to-date roads using Ordnance Survey free Open Esri Shape Data.

In particular, I want to make the map with the Contour Lines and Stepped Hypsometric Tints that John Bartholomew used - Bartholomew used non-uniform contour interval that varied with the height above sea level.

Is it possible to make Custom Contours with variable intervals in feet from DEM or other Elevation Data, or is QGIS restricted to one Contour Interval in metric units only?

The reason for this is a personal subjective one - I feel that elevation in feet gives a better 'feel' of the relatively modest range of elevation in the British Isles than contours in metres do.

If QGIS can't do it, is there another mapping application that can?


texts: Determine if ArcGISDynamicMapServiceLayer is from a Tiled Map Service in ESRI JavaScript API

texts:
For some of our JavaScript mapping applications, we load a Tiled Map Service as a ArcGISDynamicMapServiceLayer. When we use the PrintTask and want to print a map at a scale in between the tiled map scales, this works great.

For a map Table of Contents widget, however, I need a way to determine if an ArcGISDynamicMapServiceLayer has been loaded from a tiled or dynamic map service. Any recommendations?

I'm using ArcGIS JavaScript API v.3.6.


texts: Calculating statistics for multipolygon in Google Earth Engine?

texts:
I have a shapefile which contains many crop plots, around 2000. I need to calculate NDVI (normalized difference vegetation index) and then calculate statistics for each one of the 2000 plots separately.

Is there any short way to make Google Earth Engine calculate the statistics for each plot without saving each plot as a shapefile and uploading it to the system?

This is the code I have now:

my end goal is to create in the end one table which will have different statistics calculations (which I still haven't added) for each plot, ordered by the dates of the images from the Image Collection.

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2015-06-23', '2019-07-02')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
// print(clippedCol);//here I get the error messege "collection query aborted after accumulation over 5000 elements
// print(dataset,'dataset');//the same error here

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

// // Test the addNDVI function on a single image.
// var ndvi1 = withNDVI.select('NDVI').mean();

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

// // Test the addNDVI function on a single image.

// var listOfImages =(withNDVI.toList(withNDVI.size()));
// var firstImage = (listOfImages.get(0));
// var secondImage =(listOfImages.get(1));

Map.centerObject(geometry);
print('NDVI Images',withNDVI.select('NDVI'));
// Map.addLayer(firstImage,colorizedVis,'test1');
// Map.addLayer(secondImage,colorizedVis,'test2');

//Long-Term Time Series

print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.mean(), 30));


texts: Sort before/during identify duplicate values

texts:
I'm using ArcMap and I'm trying to find (and eliminate) features based on duplicate values.
I found the following code:

uniqueList = []
def isDuplicate(inValue):
  if inValue in uniqueList:
    return 1
  else:
    uniqueList.append(inValue)
    return 0

and

isDuplicate( !FIELD!)

which works great, however the problem I have is that its not always identifying the duplicates I want to eliminate. (Some are right, some are wrong.)

If there is a way to sort by another field as part of this code, that would solve the problem.

Is there something I can add to the code above to sort by another field (shape_area smallest to largest)

Obviously I've tried sorting the attribute table first, but it doesn't seem to affect the duplicate code.


texts: Rotate functions in "Rectangles, Ovals Digitizing" as well as "Advanced Digitizing" distort shape dimensions

texts:
I am trying to rotate squares via the rotate tool in the "Rectangles Ovals Digitizing" Plugin in QGIS. The dimensions/proportions of the sides, however, are distorted whenever I do so. Any idea of why this is occurring? Could it be an issue with my CPS?

The purple square would be the original square, whereas the yellow polygon is the rotated form.

Any suggestions on how I would keep the "square" dimensions of the original shape?


texts: I can not see geometry column on SQL Executor of FME

texts:
I have build a simple workflow in FME 2014 (as shown in the image) to get data from an excel file and store them in a postgres database. For this purpose I use SQL Executor.
I manage to pass all the data to the database but for some reason I can not see in the dialog box where I build the SQL Statement (in the properties of SQL Executor) the geom column which I have in my database table.
Why this happens?


texts: Creating color relief from GeoTIFF using gdaldem returns font name

texts:
I'm trying to generate a color relief map from a SNODAS greyscale geotiff file scaled from 16-bit to 8-bit.  Given the type of data and the scaling, it seems the pixel values are set at 1-5, with 5 being nearest white and 1 being nearest black.  I'm trying to create a color relief from this using gdaldem using this:

gdaldem color-relief [inputfile.tif] [color.txt] [outputfile.tif]

However, I am returned this:

ERROR 1: Unknown color : ArialMT;}

This appears to be the font used in the text file, rather than the values I've put there, which is this:

0 nodata
1 blue
2 yellow
3 orange
4 red
5 white

Any ideas here?


texts: Can ogr2ogr be used to append multiple DBF files? It works with SHP files, but I can't see how to do it with DBF files on their own. It ought to be possible, considering DBF is a subset of SHP.


texts: Problem in saving a file in a specific output folder

texts:
I have a very simple script that I'm testing. It has always been working so far, but now I'm running in an unexpected error. When I try to save my files in a folder different than the environment one, it outputs this error:

    outSnowmask.save (newpath_01+ '\\'+ 'snow_'+ MODIS)
UnicodeDecodeError: 'ascii' codec can't decode byte 0xe0 in position 47: ordinal not in range(128)

My code is:

import arcpy, os, sys, shutil

newpath_01= r'C:\\Users\\unimi\\Documents\\Umberto\\Università\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_01'
if not os.path.exists (newpath_01):
    os.makedirs (newpath_01)

arcpy.env.workspace = 'C:\\Users\\unimi\\Documents\\Umberto\\Università\\PhD\\Karakoram\\MODIS_Pakistan\\MODIS_Clip\\Sample_test'
arcpy.env.overwriteOutput= True
from arcpy.sa import *

listMODIS = arcpy.ListRasters ()
print listMODIS
arcpy.CheckOutExtension ('Spatial')

for MODIS in listMODIS:
    outSnowmask = ExtractByAttributes(MODIS, "VALUE = 200")
    outSnowmask.save (newpath_01+ '\\'+ 'snow_'+ MODIS)


texts: Reverse Geocoding within a Radius with Geotools using Natural Earth Shapefiles

texts:
The problem I'm trying to solve is: given a latitude, longitude, and a radius, get all the countries contained within in an area with that radius. My approach is to first, buffer my point by a radius using the bufferPoint method (https://gis.stackexchange.com/a/311333/176526). Next, I modified the lookup method (https://gis.stackexchange.com/a/182469/176526) work with Geometry coordinates rather than a point:
public List<SimpleFeatureCollection> lookup(Geometry g, SpatialIndexFeatureCollection countries) {
    GeometryFactory gf = new GeometryFactory();

    Coordinate[] bufferedCoordinates = g.getCoordinates();
    Point[] points = new Point[bufferedCoordinates.length];

    int i = 0;
    for (Coordinate c : bufferedCoordinates) {
        points[i] = gf.createPoint(c);
        i++;
    }

    List<SimpleFeatureCollection> simpleFeatures = new ArrayList<>();
    for (int j = 0; j < points.length; j++) {
        Filter f = ff.contains(ff.property("the_geom"), ff.literal(points[j]));
        simpleFeatures.add(countries.subCollection(f));
    }

    return simpleFeatures;
}

Finally, I iterate over the list of features:
    int i = 0;
    for (SimpleFeatureCollection simpleFeature : simpleFeaturesList) {
        SimpleFeatureIterator itr = simpleFeature.features();
        try {

            while (itr.hasNext()) {
                SimpleFeature f = itr.next();
                //System.out.println(f.getAttribute("NAME"));
                i++;
            }
        } finally {
            itr.close();
        }
    }

I do realize this is sloppy, but problem is that it's also wrong - given a large enough radius, none of the Geometry's coordinates will fall into some smaller country contained within the radius.
How do I do this correctly?


texts: How to unpolygonize and remove outermost border?

texts:
QGIS has "polygon to line" but this isn't exact reverse engineering of "line to polygon". I want no overlapped border closed lines, just segmented (on each intersected point) lines.

Additionally I want to remove to outermost boundary (if it can be automatically done).

Example: in the case of a kite shape (4 triangle polygons) will turned into inner skeleton cross lines (4 segments)


texts: Difference operation between selected features and the rest of the features of a same layer

texts:
I have got a shapefile of overlapping polygons that have an attribute "priority" (values range from 1 to 5).

I would like to do a succession of "difference" operation in order to create a new shape that don't contains any overlap based on the field "priority".

For example the polygon features with priority 1 will "erase" all the ones it overlap, the features with priority 2 will "erase" all the ones with priority 3,4,5 and so on...

My initial idea was to:

dissolve based on the "priority" field
select a feature and do a difference

I have tested the code below to see if my strategy would be okay but it doesn t work

# I have successfully dissolved my layer based on the "priority" field
layers = iface.legendInterface().layers()
layer = layers[0]

features = layer.selectedFeatures()
layer.invertSelection()
inverse = layer.selectedFeatures()

for f in inverse:
    features[0].geometry().makeDifference(f.geometry())

What is the way to achieve this?


texts: Unable to Determine Error in QGIS Log Message Panel

texts: A plugin is built for our road data network in qgis which basically splits roads at intersection and use grass v.clean algorithm, this tool runs smoothly until now, but now it's not splitting road, from snip can anyone gives me idea what is wrong, QGIS version 3.4 stable and OS is Windows 10.


texts: What is a valid WMTS Capabilites TileMatrixLimits element if the TileMatrixSet MatrixWidth or MatrixHeight is 1

texts:
In a OGC WMTS 1.0.0 capabilities document, part of the TileMatrixSetLink is the  TileMatrixSetLimits element.  It contains a list of TileMatrixLimits elements, which in turn specify the dataset extent in terms of the min/max column/row.  There must be one TileMatrixLimits for each TileMatrix in the TileMatrixSet.

The TileMatrixLimits element MaxTileRow is listed as:

Maximim tile row index valid for this layer. From minTileRow to matrixWidth-1 of the tileMatrix section of this tileMatrixSet

The Xml schema document lists the MaxTileRow type as positiveIntegers, which exclude zero.

Ignoring the fact that the description references the matrix width instead of height, how is a TileMatrix with a MatrixHeight equal to one handled?  From what I can tell, this can't be represented by a TileMatrixLimit, which means that the whole TileMatrixSetLimits element has to be removed.

Based on the description of MaxTileRow, shouldn't the Xml type be nonNegativeInteger?

Looking at the WMTS 1.0.1 schema, MinTileRow and MaxTileRow are both type PositiveInteger.

http://schemas.opengis.net/wmts/


texts: Creating “oblique bounding box” with maximum width/height ratio in QGIS?

texts:
In QGIS I would like to create bounding boxes that align with the polygon footprint.

When using "Oriented minimum bounding box" it skews the bounding possibly based on the world aligned positions.

When using "Bounding Boxes" it has a similar issue using a true min max position.

The desired result I am after would be this:

There is a similar question already posted here. It offers a script that seems to have a good answer. I would like to know how to achieve the same thing in QGIS. I have tried to find an equivalent of ST_Rotate etc but have not found anything specific.


texts: hardcoding port number in Arcgis License Manager

texts:
In the Arc10 Q&A for the 2011 user conference it says a performance gain may be realized by hardcoding the license manager to port 27000, and then instructing the firewall and Windows Defender to ignore that port.

That same doc references the online help for the specifics of how to implement that on the license manager server, namely edit services.txt to look like the below (the trailing numbers on first two lines used to blank):

SERVER this_host ANY 27000
VENDOR ARCGIS PORT=1234
FEATURE ACT ARCGIS 1 permanent 1 7ED49106D630

This is all well and good to this point, but I'm confused about what is happening with the server port and the vendor ports. Which way does the communication flow go? If the ArcGIS Vendor port is hardcoded, then does every client need to be updated to explicitly match that port? (e.g. do I have to visit every client machine and set their ArcGIS Desktop Adminstrator to use port 1234? or is it 27000?) Do I set the firewall to ignore one port or both ports? in or out? And is #1234 actually a safe number, or range of numbers?


texts: How to set .gdb as Python script tool parameter?

texts:
I was thinking the Geodataset (as shown in figure 1) would be the right choice if I want to use .gdb as a parameter when using script tool, but it didn't work since I couldn't add my gdb to the box, shown in figure 2. What is the right parameter type I should use here?

Using 10.3.1 ArcGIS Advance License.


texts: What does number of shapes does not match the number of table records mean when trying to repair a data source?

texts: I'm trying to repair a data source on arc from a .shp file. Every time I match it with the correct .shp file it display an error message stating "Number of shapes does not match the number of record tables". It also will not let me add the desired .shp file to the map either. Is there a way to repair this connection


texts:
I have a layer (base layer) with MultiLineString geometries and now I need to project the layers (overlays) in the base layer from PostGIS using OpenLayers.

However, I got an idea how to project table data source into OpenLayers. I am trying it on topology schema and I am having my layer data in edge_data, node, faces, relations tables. Can anyone guide me how to project my topogeometry into OpenLayers?

And moreover, I am not using any servers except PostgreSQL.


texts: Why do some shapefiles (e.g., state-level in US) have different default "zoom levels" to others?

texts:
Why is it that this shapefile of Washington state, when loaded into QGIS or via D3, appears large and zoomed in, like below:

Meanwhile, others, such as this Zillow neighborhoods shapefile ( also for Washington state), when loaded, looks thus:

While the first shapefile was loaded at the correct, state scale, the second, neighbourhoods file was loaded at an incorrect scale — a national one. For reference, here's what it looks like at the appropriate scale:

Is there any way to control this behaviour? Can I ensure that the Zillow shapefile of Washington state's neighborhoods loads at a larger, more zoomed in scale?

edit: this happens with geojson and topojson too.


texts: Cesium camera.setView to rectangle clipping the rectangle

texts:
My application allows users to add GeoJSON objects to a cesium map.  If the GeoJSON has a bbox, it will set the view to the bounding box:

const destination = Cesium.Rectangle.fromDegrees(
    bbox[0], bbox[1], bbox[3], bbox[4]);
this.viewer.camera.setView({destination});

My naïve expectation is that the entirety of the GeoJSON object would be shown in the view, but that is not what I am experiencing.  It appears that if the bbox has the same (or quite similar) aspect ratio to the div containing the cesium map then it behaves as I want, but if not then the view only shows a portion of the GeoJSON object.  For example, this is what I am seeing:

And this is what I would expect to see:

Note that the first image is exactly as rendered after setView, the second image is shown after manually zooming out, and the blue/white stripped box is the bbox drawn to demonstrate the issue.

Further investigation seems to indicate that if the bbox is width constrained, then Cesium will zoom out to include the entire bbox, but not if height constrained.  To say it another way, if the bbox is tall and skinny, I have a problem, but if short and fat, I don't...


texts: Install GDAL python package in virtualenv on mac

texts:
I'm trying to instal the GDAL python package, preferably in a virtual environment, but when I try "pip install gdal", I get the following error:

extensions/gdal_wrap.cpp:3085:10: fatal error: 'cpl_port.h' file not found
    #include "cpl_port.h"

I have GDAL version 1.11.4 installed on my Mac. I installed it via Kyngchaos, and it is located at /Library/Frameworks/GDAL.framework/Versions/1.11/Resources/gdal

the cpl_port.h file is located in /Library/Frameworks/GDAL.framework/Versions/1.11/Headers

I've read several related topics but they're all on ubuntu and using older versions of pip that have commands like "--no-install" etc.


texts: Unable to create ArcGIS network dataset perhaps due to size?

texts: I have a large (several hundred thousands vectors, 800-900) from which I am trying to create a network (road) and am unable .  Don't sure if its an issue with ArcGIS or limitations with my computer.  That said are there any alternative methods of creating LARGE network data sets.


texts: Grid lines and text together with Atlas. Problems with some scales for each map

texts:
I'm having trouble using the "Grids" option in the QGIS 6.3 composer layout.
If I define the Interval units in Centimeters the coordinates on the map (text) are NOT referenced to the coordinate systems, in my case it is 4674, they are relative to the intervals and not in real coordinates.

If I change to MAP UNIT, the coordinate texts conform to the 4674 system, but in combination with the use of the ATLAS option, which automates several maps with different scales, I get a GRID variation on the map, which on a certain scale(1/500) is with good spacing and for smaller scales (1/3000) is very close spacing, overlapping the coordinates lines and text.

What could it be to do to correct? It would be possible to maintain grid spacing in centimeters, but would the grid text be referenced to the coordinates system? Or, allow the use of formula so that the grid spacing varies according to the map scale controlled by ATLAS?
Any suggestion?


texts:
I've downloaded some .asc files from the UK Government. They are 1m resolution LIDAR tiles. Apparently.

When I load them into QGIS, they appear with the nice thumbnails and everything. But I get errors about no CRS and defaulting to CRS EPSG:4326 WGS84. I think that is the standard model for the globe?

When I try to load the .gpx file for my route, it is no where to be seen. I think the .asc files dont have coordinates or something? I was expecting the trace to appear over the top so I could pull out surface elevations.

Any ideas?


texts: Converting raster to integer values and exporting as 16 bit using QGIS

texts:
I am trying to export a DEM raster as a 16 bit image using the Raster>Conversion>Translate drop down menus in QGIS.  I have seen somewhere I need to save it with a -ot tag to turn it into a 16 bit image but I cannot work out how to do that.  If I type in -ot into the end of the output file name it does not work . I know I am not getting it.
Alternatively, how do I do it via command line?


texts: Setting default raster statistics extent in QGIS?

texts:
I recently found the setting in QGIS to rescale raster symbology using the current extent, which is really handy (Layer styling > Min/Max Value Settings).

My question is - can the default extent be changed? At the moment, when you add a new raster to the map, the statistics extent is set to Whole raster. I would like the default to be Current canvas.


texts: Is there a source for historical shapefiles of Europe in the modern era?

texts:
Is there a source for historical shapefiles of Europe in the modern era, i.e. starting ca. 1400?

I know that "where can I find historic border data" questions were asked multiple times but I couldn't find an answer to my specific question.


texts: Labeling Layers in Data Driven Pages chops those outside of page?

texts:
I am trying to set up a soil layer to be clipped to the extent of my current data driven page.

In Data Frame Properties I have set to clip layers to current data driven page extent.

The flaw with this is that labels for this layer that are outside of the page get chopped off rather than adjusting to the new extent.

i.e. a soil that extends outside of the DDP and is labeled on the outside doesn't show up on my map.

Does anyone know how to change the label properties so that these labels come into the extent of the clip?


texts: inserting points in postgis with pdo (php)

texts:
I have a problem with sql insert statement using PDO (php, postgis, postgres). I have tried to use different combinations, but none of them work. Here is what I ve got:

$latitude= 46.07175;
$longitude=16.866976;
...
$sql = "insert into table (geom, name, type) values ('(st_setsrid(st_makepoint('|| $longitude ||',' || $latitude || '),4326))'::geom, :name, :type");


texts: How to use GDAL.Grid() with vrt file in C++/ C#?

texts:
Referring to the GDAL_GRID C++ documentation here.

gdal_grid [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
          CInt16/CInt32/CFloat32/CFloat64}]
          [-of format] [-co "NAME=VALUE"]
          [-zfield field_name] [-z_increase increase_value] [-z_multiply multiply_value]
          [-a_srs srs_def] [-spat xmin ymin xmax ymax]
          [-clipsrc <xmin ymin xmax ymax>|WKT|datasource|spat_extent]
          [-clipsrcsql sql_statement] [-clipsrclayer layer]
          [-clipsrcwhere expression]
          [-l layername]* [-where expression] [-sql select_statement]
          [-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]
          [-a algorithm[:parameter1=value1]*] [-q]
          <src_datasource> <dst_filename>

src_datasource:
Any OGR supported readable datasource.

It says that one needs to pass in a GDALdataset as an input.

However, for the example given on the same page, it is passing in a file ( with the *.vrt) as input.

gdal_grid -a invdist:power=2.0:smoothing=1.0 -txe 85000 89000 -tye 894000 890000 -outsize 400 400 -of GTiff -ot Float64 -l dem dem.vrt dem.tiff

Obviously the example is geared towards the command line utility, not the actual API usage.

I want to use the API and not the command line utility for GDAL_GRID, so how can I call the gdal_grid API properly, with .vrt dataset? Is there any such example in C++ or C#?

PS:

There are a lot of examples already in python on how to pass in a .vrt file as input, but this is not what I'm looking for. Just to clarify.


texts: Problem displaying Line Pattern Fill style on Qgis2web

texts:
I need to export a categorized map of a region from qgis via qgis2web into an html page. The first test went well, but the problem is that I can't visualize this style applied to some municipalities on Qgis2web.
This is the result on QGIS:

And this is what I see on Qgis2Web:

To create this style I used both Line Pattern Fill and Simple Fill within the categorized symbology
Is it possible to fix and display it in the preview and exported file ?


texts: How to batch open and export to PDF folder of map packages (*.mpk)?

texts:
I have a directory full of ArcGIS map package files (mpk) that I need to export to pdf. Is it possible to automate it so I don't have to open and export each one manually?

I need to open each map package, switch to layout view, export to pdf, close it, and repeat with the next map package file.


texts: Unable to draw a polygon due to "The feature cannot be added..."

texts:
The feature cannot be added because it's geometry collapsed due to
  intersection avoidance.

see also http://i.imgur.com/xiiI7NE.png

What happended: Let's say I have 2 layers. To the first I enabled the snapping tool with avoid Intersections. So far so good.

Now I added a polygon to the layer 2. That layer has nothing enabled in the Snapping tool. But as you can see from the error above, it's broken. So I'm not able to create a polygon on top of the layer 1. Is this intended? I mean where does the avoid intersection option make sense if this is intended?

This is more a general question, because with disabling the avoid intersection I was able to create a polygon on layer 2.


texts: Index accelerating complex feature queries and geometric computations?

texts:
As far as I know, the index on geospatial database only optimizes lookup of features, not computation themselves, eg, it will acceleration intersection between a large number of polygons with a small number of vertices, but will not accelerate much with a lower number of polygons with the same amount of total vertices.
Besides splitting the polygons, am I missing something on geospatial databases practices to achieve great performance on operations (eg, intersection, nearest) against polygons with a large number of vertices ?


texts: Shapefile to KML not rendering in Google Earth

texts:
I have this shapefile which does not open in Google Earth when converted to kml.

GCS: WGS1984

PCS: UTM, zone 41N

I have been trying it with GCS only but it still does not work. I have been using the define projection and Feature--> project tools but it did not work either... other kml files, on my system, are working.

Can someone indicate where I'm getting it all wrong. I have seen and tried solutions given in similar questions:

How to convert coordinates from meters to decimal degrees?

but they don't help, but, rather, I lost the file and see an empty canvas. Google Earth positions it somewhere in the Arctic but the polygon etc are not visible.

The coordinate system is in meters which I normally have in decimal degrees.


texts:
I have a polyline, and i want to add point (vertex) on a specific distance, like two meters for example.

How can i do that?


texts:
I try to make a prediction in raster using linear regression based on past rasters but I can not capture the lm() model in order to feed it to predict().

#
# calculate regression in a raster stack (time series)
# and predict.
#
require(raster)
require(rgdal)

# make the list of rasters
rasters <- list.files(pattern='*.tif$')

# make the raster stack
s <- stack(rasters)

# crop the stack to the extent of Sicily
sicily_ext <- c(12, 16, 36.5, 38.5)
sicily <- crop(s, sicily_ext)

# make a time variable (to be used in regression)
time <- 1:nlayers(s)

# run the regression
fun <- function(x) { if (is.na(x[1])){ NA } else {lm(x ~ time)$coefficients[2] }}
x2 <- calc(sicily, fun)

#predict to a raster
predicted <- predict(sicily, x2, progress='text')

# END

The exact output is:

Loading required package: raster
Loading required package: sp
Loading required package: rgdal
rgdal: version: 1.0-6, (SVN revision 555)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
 Path to GDAL shared files: /usr/share/gdal/1.11
 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.1-1

0%Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "data.frame"

and the resulting x2 raster is shown here


texts: QGIS graphical modeller and statistics tables

texts:
I have been able to use a processing tool available in QGIS, 'Statistics by Categories'to obtain some statistics for lengths of pipes of different materials within a polygon area. I am only actually interested in the sum of the lengths of pipes within each of the material categories.
I had therefore hoped to be able to either select the statistical measure of interest and ignore others within the tool while setting up the process, or I have tried to use the 'drop fields'algorithm on the result to generate a further result that only has my measure of interest included (but this only works on vector layers). Neither method does what I am hoping for.
Would anyone know of an appropriate processing tool or technique to achieve this?


texts: How to clean the postgres instance while setting cartoDB locally in ubuntu

texts:
I am following https://github.com/danseely/cartodb-install/blob/master/DEV-INSTALLATION.md for setting cartoDB locally in ubuntu12.04.

While installing cartoDB-postgresql, the second install check of cartoDB-postgresql fails in test_ddl_triggers and also in CDB_quota_test.

It is mentioned to change the shared_preload_libraries = 'schema_triggers.so' in postgresql.conf if test_ddl_triggers and restart the postgres by service postgresql restart but still the error remains same after making changes and restarting.

It is also given to clean postgres instances in https://github.com/cartodb/cartodb-postgresql/ to test installation.

How to clean the instances and test the installation for error free?


texts: If I were to draw a 1km square on the earth what would the angles of the corners be?  Obviously they would be a little greater than 90-deg but how would I calculate the angle exactly?


texts: How to make a custom Arc toolbar part of the default document?

texts:
At ArcMap 10.7.1, I have created this custom toolbar, which contains a single tool:

I need that toolbar to "stick", so that it always appears when I open a new or existing map document (I think that means my new toolbar becomes part of normal.mxt?).

However, this toolbar does not appear when I subsequently open a map document.  According to this ArcGIS help document, I need to have Customize > Customize Mode > Options > "Save all customizations to the document" checked off, as shown here:

However, I am unable to get it to "stick", no matter what combination of boxes I check and uncheck; the custom toolbar does not appear when I open a map document.

This issue was brought up in SE GIS over three years ago in 10.1.  Here we are six-plus releases later... is it still true that you cannot make a custom toolbar part of the ArcMap default view???

EDIT:

Here's a quote from the 10.7 Help page (my emphasis):

When you create a custom toolbar—in ArcMap, for example—it will save with the current map document by default. If you want to save a toolbar in the Normal template so that all documents see it, uncheck the Create new toolbars and menus in document check box on the Options tab on the Customize dialog box.


texts: How to restrict level crossings in route?

texts:
I want to evaluate if it is useful to bypass a level crossing through a nearby tunnel or overpass when the LC is closed or if it is worth waiting.
I want to use the OSM data for realistic scenarios and calculate the most time efficient route for closed and opened LC.
Is it possible to solve this problem with pgRouting?
Or can you give me some hints where to start?


texts: Saving Metadata to DB (QGIS 3.4 or higher)

texts:
Is it possible to save metadata of a QGIS 3 layer in a database (PostGIS), similar to the function in QGIS 2.18 "save layer style to database".
I don't want to save it in the project or as a file.


texts: Issue displaying UTM coordinates in QGIS 3.6.1

texts:
I am having trouble properly displaying my .shp layers in UTM (QGIS 3.6.1). Whatever I do, the coordinates always display as lat/ long, and measurements are way out. I have re-saved all my .shp layers using the CRS GDA94 MGA zone 54, as well as setting my project CRS to GDA94 MGA zone 54, too.

Am I missing something?


texts: How to replace QGIS featureAtId() with getfeatures()?

texts:
I'm trying to upgrade a script made for QGIS 1.8 to work at version 2.14

One of the problems I face is that the featureAtId function at 2.14 is now obsolete. Original codes looks like this.

        for feature in coastal_10k_layer.dataProvider().fields():
            coastal_10k_layer.featureAtId(feature, f)
            geo = f.geometry()
            if geo:
                geometry = QgsGeometry(geo)
                self.geometries_10k.append(geometry)

I tried something like this as official site recommended and doesn't work.

        for f in coastal_10k_layer.getFeatures():
            geo = f.geometry()
            if geo:
                geometry = QgsGeometry(geo)
                self.geometries_10k.append(geometry)

Is there any way to do this?


texts: Retrieving Landsat cloud percentage without downloading the data

texts: Is there any way to get the cloud percentage information from metadata of Landsat data without downloading it? I would like a Python solution if one is available.


texts: I have a series of XY coordinates. I want to use these points to draw a line in ARCMAP using arcobject10.1 using c#. I have saved this in a table inside geodatabase.Am i supposed to create a feature class for this purpose or should i use graphics. I would appreciate a hint in this direction.  I have been working all day without any success so help is needed. Please note as i mentioned  i am looking for arcobject solution not python, :)


texts: Set Root Mean Square Error (RMSE) as feature properties

texts:
Using Landtrend library in Earth Engine, I can use a property that comes as default to calculate the square mean error (RMSE) and, then use this property to apply the current formula for this model:

  (Math.round(ltData.rmse * 100) / 100).toString()

I just need to select 'rmse' property inside the 'LTRESULT'. First, I need to reduce to the level of each one of the fifth samples (samplesVegPrimaCinquenta).

LTresult.select("rmse").reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: samplesVegPrimaCinquenta,
  scale: 30,
}).evaluate(function(ltData){ ...

I want to return this value (RMSE) as a property inside each feature in the FeatureCollection, but I can't do this.

var triplets = collection.map(
  function(image) {
  return image.reduceRegions({
    collection: samplesVegPrimaCinquenta,
    reducer: ee.Reducer.mean().setOutputs(image.bandNames()),
    scale: 30,
  }).map(
    function(feature) {
    return feature.set({
      'imageID': image.id(),
      'date': ee.Date(image.get('system:time_start')).format('Y-MM-dd'),
       // Calculate rmse and set as propertie inside each feature in featureCollection
      'Fit RMSE: ': (Math.round(ltData.rmse * 100) / 100).toString()
    })
  });
}).flatten();
print(triplets.limit(5),'triplets')

Here's the LINK code, the script is a little big so, this part it is in line 178-208


texts: Parsing list object file paths for file name

texts:
I am trying to automate a naming process that takes the file name from a list of file paths to name other folders, geodatabases, zipfiles, etc.

The problem I'm running into is that, I assume, I am trying to run string functions on list object and I can't figure out how to get the script to do what I need it to do.

For example, I have the file path C:\TEMP\Durham.shp. I want to name a geodatabase, folder, and zipfile "Durham". Instead I get a folder and geodatabase named "C".

Can someone take a look at my script and tel me what I'm doing wrong?

The following script is that last incarnation of several attempts.

import arcpy, os
import zipfile
from arcpy import env
from os.path import basename

BASELAYERS = arcpy.GetParameterAsText(0)
CLIPLAYERS = arcpy.GetParameterAsText(1)
OUTFOLDERMASTER = arcpy.GetParameterAsText(2)
FILETYPES = ["*.shp", "*.dbf", "*.shx", "*.prj", "*.sbn", "*.sbx", "*.xml"]

env.workspace = OUTFOLDERMASTER
for LAYERS in CLIPLAYERS:
    LAYERSSTR = ''.join(LAYERS)
    BASENAME = LAYERSSTR.split("\\")[-1]
    FOLDER = os.path.join(OUTFOLDERMASTER, BASENAME)
    os.makedirs(FOLDER)
    PGDBNAME = "NCFlood_Effective_" + BASENAME + "_PGDB.mdb"
    arcpy.CreatePersonalGDB_management(FOLDER, PGDBNAME, "9.3")
    PGDB = os.path.join(FOLDER,PGDBNAME)
    for FILES in BASELAYERS:
        FEATURENAMESTR = ''.join(FILES)
        FEATURENAME = os.path.basename(FEATURENAMESTR)
        SHAPEFILENAME = FEATURENAME + ".shp"
        SHAPEFILE = os.path.join(FOLDER, SHAPEFILENAME)
        arcpy.Clip_analysis(FILES, LAYERS, FEATURENAME)
        arcpy.Clip_analysis(FILES, LAYERS, SHAPEFILE)
    ZIPFILELIST = []
    for TYPE in FILETYPES:
        ZIPFILELIST.extend(glob.glob(FOLDER))
    ZIPFILENAME = os.path.join(FOLDER, "NCFlood_Effective_" + FOLDERNAME + "_SHP.zip")
    ARCHIVE = zipfile.ZipFile(ZIPFILENAME, "w")
    for ZIP in ZIPFILELIST:
        ARCHIVE.write(ZIP, os.path.basename(ZIP), zipfile.ZIP_DEFLATED)
    ARCHIVE.close()

I have also tried

    BASENAME = os.path.basename(LAYERS)
    FOLDERNAME = os.path.splitext(BASENAME)[0]


texts: Applying sorted sequential numbers to only certain features using ArcPy

texts:
I would like to take sequential numbers up a level and only apply it to features of a certain type. Say I have two feature types, black and white features, and it's denoted by the field "Type". I want my sequential numbers code to skip all white features and pick up where it left off with the last black feature. I'm sure it's some type of If/Then statement but i can't figure it out. Here's the basic code where it returns sequential numbers by my sort field, which works just fine:

import arcpy

sortFeat = r"P:\path\Receivers.shp" #featurepath
sortField = 'Ycoord' #Base Field to sort, which is
idField   = 'R_ID' #Field to populate sequential numbers
rec=0

def autoIncrement():
    global rec
    pStart = 36 #number of features
    pInterval = 1
    if (rec == 0):
        rec = pStart
    else:
        rec -= pInterval
    return "R" + str(rec)

rows = arcpy.UpdateCursor(sortFeat, "", "", "", sortField)

for row in rows:
    row.setValue(idField, autoIncrement())
    rows.updateRow(row)
del row, rows

print "Finish"

So how can i make this only apply to certain attribute type and not all within the layer?

I am using Python 3, but this currently only works in a 2.7 environment.


texts: Feature splitting when pasted into another shapefile

texts:
Using QGIS I was trying to copy a 2 part feature into another shapefile as a single entity. I was able to create a single feature

the way I wanted it. I then attempted to copy and paste it into the ultimate destination layer but it then split/exploded into several separate features.

as can be seen by each one being labelled. The individual features are all separate in the attribute table. The other feature in the screen shot is as I expected to see the new feature.
I have repeatedly attempted various approaches but when the result is the same. I am sure I must have done this successfully before.
I have only recently started using 3.4.15-Madeira on Mac Catalina so the problem might lie there.


texts: Calculating mean upslope aspect from each cell in DEM using ArcPy?

texts:
I am trying to write a python script to calculate the mean aspect upslope from each cell in a DEM.  I have the general workflow down, but my final output is restricted to angles between 0 and 90 degrees, which is not correct for the input data.  The general procedure to average an angle is:

Calculate the sine and cosine of each angle.
Sum the sines and cosines.
Use atan2 on the sums to find the mean angle.

I have successfully used this technique to find focal and zonal mean aspects, and am trying to adapt my code to calculate upslope aspect.  I am summing up the sine and cosine rasters using a weighted flow accumulation, which is the only difference between this and my zonal/focal workflows.  I have included the script below.

Does anyone have any idea what's different about flow accumulation that causes this, or is there something else in my code that I'm missing?

arcpy.CheckOutExtension("spatial")
import math
from arcpy.sa import *

#collect input parameters
inDEM = arcpy.GetParameter(0)
flowDir = arcpy.GetParameter(1)
outMean = arcpy.GetParameterAsText(2)
scratch = arcpy.GetParameter(3)

#set workspace and snap raster
arcpy.env.workspace=scratch
arcpy.env.snapRaster = inDEM

#Calculate aspect and set flat cells to NoData
rawAspect = Aspect(inDEM)
nullFlat = SetNull(rawAspect,rawAspect,"Value < 0")

#convert aspect to radians and calulate cos/Sin
Radians = Times(nullFlat,0.01745329)
cosAsp = Cos(Radians)
sinAsp = Sin(Radians)

#sum upslope Cos/Sin rasters, use ATan2 to average
cosAccum=FlowAccumulation(flowDir,cosAsp)
sinAccum=FlowAccumulation(flowDir,sinAsp)
ArcTan = ATan2(sinAccum,cosAccum)

#convert mean aspect back to degrees and save output
meanAspect = Mod(360+ArcTan*(180/math.pi),360)
meanAspect.save(outMean)


texts: Boolean If then statements in ModelBuilder and Repeating model?

texts:
I have a large, several step model that iterates through feature classes in a folder taking several time-stamped lat long points and populates several fields in  the attribute table.  After the model is complete, I am looking to establish and if/then statement based on the derived 'KMPerHour'.  If any records are above 40KPH, I need to remove the points that are above the that threshold, and repeat the entire model using the subset of the data.  Once all segments are below the 40 kph threshold, the model will stop.  I know that I will need to use a while loop, but am unsure how to do so.  I am attaching two screenshots.  One is a screenshot of the end of the model, and the other is the beginning of the model with the iteration.  if someone could help with the next steps to establish the if then statement:

If all records of 'KmperHour'<40 then stop the model
else delete all points above 40 and use the new dataset as the input of the model.


texts: QGIS+PostGIS set atlas extents to include atlas feature and point

texts:
I have a QGIS print layout with an atlas layer for which I would like the map extents to show the current atlas feature and a nearby point. The layer provider is postgis.

It seems the only way to do this is provide separately the xmin,xmax,ymin,ymax via expression to qgis. I've also discovered it behaves oddly if the aspect ratio of the bounding box is different than the map item so I need to account for that too.

What I've tried so far:

Get a bounding box from PostGIS scaled by 10% height

SELECT t.territory,
 st_expand(t.bbox, ((st_ymax((t.bbox)) - st_ymin((t.bbox))) * (0.1))) AS geom
FROM
 (SELECT
   parcellayer.id,
   st_envelope(st_union(parcellayer.geom,
    (SELECT pointlayer.geom FROM pointlayer LIMIT 1))
   ) AS bbox
 FROM parcellayer) t;

this gets me a box but I need to make sure it's the same aspect ratio (2:1)
I can do a bunch of complicated case statements, but is there an easier way to go about this?


texts: Creating evenly spaced grid for QGIS composer to show overlap in Map?

texts:
I would like to create an evenly spaced, overlapping grid so that I can use the QGIS ATLAS to display a large mapping extent with a little bit of overlap between each image of the Composer.

I am working in QGIS 2.18.14 and have not been able to figure this one out.

I have tried using the 'Grids for Atlas' Plug-in but it does not give the option to set a uniform Extent when using its dynamic grid option.


texts:
How can I add my 'slippy map' (the XYZ tiles) designed in carto(db) into QGIS?

I am aware of the CARTO plugin for QGIS but it cannot import a carto layer with the styling done in CARTO; it only imports the data and I would have to restyle the data again in QGIS.

The Quickservicesmap plugin for qgis only includes the CARTO base maps designed by Carto; I am trying to add a custom map that I already heavily styled in carto.

I've done some digging From Carto's maps API and think I may be able to obtain the tiles from carto. From reading that link above, I'm not sure if I have to initialize an API call to create the XYZ tile layer first and if I do, how I'd do that in QGIS (or if there's another way to do it).


texts: Applying OSC changesets with negative node/way ids

texts:
I'm trying to apply a list of .OSC changesets over an .OSM file.

The problem I have is that the .OSC files contain negative ids for nodes and ways.
I understand that the negative ids are merely placeholders and that they are only unique per OSC file/entity type.

The fact that I want to merge such chagesets creates some (expected) issues:

> osmconvert initial.osm changes1.osc changes1.osc > updated12.osm
osmconvert Warning: wrong order at node -1 in file changes1.osc
osmconvert Warning: next object is node -2

osmconvert and osmosis complain about ids not being sorted in ascending order (in OSCs the ids are -1, -2 and so on)
any item from one changeset will overwrite elements from the previous changesets that happen to have the same id

I looked at the options that osmconvert and osmosis offer for dealing with negative ids, but found nothing - it appears that they assume all files would contain real ids.

I could write a preprocessing tool that assigns positive ids, but was wondering if there's a way to deal with this situation using the existing tool set.


texts: Carto PackageManager return null on Xamarin.iOS

texts:
I have such problem: On Xamarin.IOS Package Manager returns null after invoking GetLocalPackage(packageId) in several places of the application;

What do I know about the issue:

1)  This happens immediately after downloading of map package. In other words, in other places of app this work normally;

2)  Maybe invoked not on main thread;

3)  Reproducible only on 5 and 5C and 5SE iPhones, but on 5C, 5SE in one case less than on 5;

I am supposing that it has some relation to the device performance.

Any ideas whould be usefull.


texts:
I installed PostGIS on a server and now I want to query the OSM ids referring to a name, in order to get a list of all geometries according to this name.

Is there an efficient way to query this ids in PostGIS or do I have to use a tool like nominatim?

I am new to OSM and PostGIS.


texts: Python: Extract median by polygon from multiple raster files

texts:
What I'm trying to do:

Extract median (or 50th quantile) by a single polygon from 558 GeoTIFF rasters stored in a folder

What I got so far:
Using Python from within Qgis and Qgs.ZonalStatistics I've managed to extract statistics by polygon from *.tiff files in a directory, and add each value to a new column in the above polygon shapefile (column name based on raster filename). The code looks like this:

import glob, os, qgis.analysis
vectorlayer = qgis.utils.iface.mapCanvas().currentLayer()
rasterfolder = 'D:/effect_test/'
os.chdir(rasterfolder)
for lyr in glob.glob("*.tiff"):
    qgis.analysis.QgsZonalStatistics(vectorlayer, lyr, attributePrefix=lyr, rasterBand=1).calculateStatistics(None)

What's left to resolve:

With the above code I'm getting statistics for average, rejects, max
and min while I need to extract the median (which is the same as
the 50th quantile). I've been trying to add QgsZonalStatistics.Mean
suggested in this question and also QgsZonalStatistics::Mean
suggested in the QgsZonalStatistics documentation. But it fails
telling me 'QgsZonalStatistics' is not defined. How to re-write above code to get median or the 50th quantile instead?


texts:
I have asked this question in StackOverflow, but got no useful answer - So perhaps I could get some help here.

I would like to generate a 3d Terrain using THREE.js, and I would like to overlay OpenLayers upon it. Does anyone know how to do it?

EDIT:
I tried this (based on this)

var sphere2 = new THREE.Mesh(new THREE.SphereGeometry(60, 100, 100),material2);
scene.add(sphere2);
sphere2.position = new THREE.Vector3(100,230,140);
var material3= new THREE.MeshPhongMaterial();
material3.map = THREE.ImageUtils.loadTexture('http://vmap0.tiles.osgeo.org/wms/vmap0?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=basic&CRS=epsg:4326&BBOX=-90,-180,90,180&WIDTH=500&HEIGHT=500&FORMAT=image/png');

This should generate a sphere, and map the texture as returned by the WMS call to the sphere, making a full earth - but it does not load anything - the screen remains black.


texts:
Currently i have created the layer in Geoserver by configuring a native SQL statement, there i am getting the route geometry details from DB by following query

   SELECT ST_MakeLine(route.geom) FROM (
    SELECT geom FROM pgr_fromAtoB('ways', %x1%, %y1%, %x2%, %y2%
  ) ORDER BY seq) AS route

this will give me linestring but along with this i want some other details i.e name of the node etc for that i just modified the query as

 SELECT ST_MakeLine(route.geom),text(route.name) FROM (
    SELECT geom,name FROM pgr_fromAtoB('ways', %x1%, %y1%, %x2%, %y2%
  ) ORDER BY seq) AS route

name above is nothing but the name of the route and text() will returns me the name if each route.But this gives syntax error.How to overcome this issue?


texts: Clipping is not including all parts of the input layer

texts:
I'm having trouble with clipping a certain layer in QGIS 3.10.13. As you can see in the 4th picture, not all of the layer is being included when clipped - only one part of the layer which is highlighted in yellow.
How can I get all of the input layer to be included in the clip?
I added a picture of the land-use_pg layer in the layers panel. The result of the clip is the yellow on top of the grey buffers, so it seems like when I'm clipping it is only including the 'agriculture' part of the layer instead of all of them.
I added a picture of the output of the validity check of layer 'land-use_pg'. This seems to explain why only 'agriculture' is included when clipping layer 'land-use_pg', but I'm not sure how to fix this.
I used the fix geometries tool and clipping works fine now.


texts: How to restrict buffer extent within land area only?

texts:
I need to create buffer for a point feature, but the buffer area will be restricted within the land area only.

For example consider the following:
A buffer of say 500 km radius is to be generated around Washington or Rome. However the generated buffer will be covering the land area of US(west-ward) or Europe(towards east or north) only, respectively and not the water areas.

How can I achieve this, preferably in PostGIS. The input tables will be land area boundary polygons and a point layer.


texts:
I would like to superimpose openstreetmap data over a raster image, using mapnik.
My problem is that this image is not in a projection being defined by a SRS, I have some non-linear algorithm to run to be able to match a pixel and a geographic position.

Is it possible to use some kind of custom converter for positions in mapnik ?

Another possibility would be to modify every longitude/latitude in the database to set them to the corresponding pixels, but I'm afraid that it won't play well with the internal reprojections in mapnik.

I still have the possibility to do everything by hand using the OSM API, but the result would be ugly.

The processing will be done in python.


texts:
I am using QGIS 3.16.3.
I wrote this script to loop through zipcodes and layers to create a layer and export as an image. However, I can't get it to add the map layer.
import time
from qgis.PyQt import QtGui

project = QgsProject.instance()
manager = project.layoutManager()
layout_list = manager.printLayouts()
for layout in layout_list:
    manager.removeLayout(layout)

for x in ["92007","92037","92011"]:
#    for type in ["Frontage Map", "Housing Stock Map", "FAR Map"]:
    zip = str(x)
    type = "FAR Map"

    canvas = qgis.utils.iface.mapCanvas()

    baselayer = QgsProject.instance().mapLayersByName("Positron (retina)")[0]
    layer = QgsProject.instance().mapLayersByName(type)[0]
    canvas.setLayers([layer, baselayer])
    canvas.refresh()
    time.sleep(1)

    layer.setSubsetString("SITE_ZIP="+str(x))
    layer.selectByExpression("SITE_ZIP="+str(x))
    select_extent = layer.extent()
    time.sleep(10)
    canvas.refresh()

    time.sleep(1)
    canvas.zoomToSelected()
#        time.sleep(2)
#        canvas.refresh()
    time.sleep(1)
    layer.removeSelection()
    time.sleep(1)
#        canvas.refresh()
#        time.sleep(2)

    fp = "somefilepath"

    # Generate Layout
    layout = QgsPrintLayout(project)
#        layout.initializeDefaults()
    project.layoutManager().addLayout(layout)
    layout.setName(type+" - "+zip)

    map = QgsLayoutItemMap(layout)
    map.setRect(20,20,20,20)

    #Map Settings
    map.setBackgroundColor(QColor(255,255,255,0))
    layout.initializeDefaults()
    ms = QgsMapSettings()
    map.setLayers([baselayer, layer])
##    ms.setLayers([baselayer, layer])
    canvas.refresh()

    #time.sleep(5)

    rect = QgsRectangle(layer.extent())
    rect.scale(1.0)
    layout.addLayoutItem(map)
    ms.setExtent(select_extent)
    map.setExtent(select_extent)
    map.zoomToExtent(rect)

    map.attemptMove(QgsLayoutPoint(0,0, QgsUnitTypes.LayoutMillimeters))
    map.attemptResize(QgsLayoutSize(295,210, QgsUnitTypes.LayoutMillimeters))

    exporter = QgsLayoutExporter(layout)
    exporter.exportToImage(fp, QgsLayoutExporter.ImageExportSettings())

Edit :
I've done a bit more testing and have gotten my layers to display by clicking the following settings from the Layout Item Properties Menu. (see png attached)

Set map extent to match main canvas
Set map scale to match main canvas scale
Lock layers
Lock styles for layers
Update map


texts: Speeding up Google Geocoding API Slow performance

texts:
I am using Google geocoding API service with their Java client, and my dev Key,
The problem is when I run up to 40 requests it does slow to 600ms. I am using it in a webservice that is supposed to run faster than that. How can I speed things?
I already set up the QPS in my code but it doesn't improve performance.

GeoApiContext context = new GeoApiContext().setApiKey("AI....");
context.setQueryRateLimit(100);
GeocodingResult[] results =  GeocodingApi.geocode(context,Adresse).await();

What do you suggest to improve the performance?


texts: Using QGIS to workaround arcAdvanced license, one sided bufffer

texts:
I need to export data from a file gdb to a shapefile or PostGIS format, run a script that will create one-sided buffers in QGIS for the sides of streets whenever a point feature is found on that street, and then import the resultant polygons back to the original file gdb. Can this be reasonably done with pyqgis calling the Single Sided Buffers (and offset lines) for lines tool in the GDAL/OGR geoprocessing library? My sample code thus far is just trying to import the qgis module but my citrix vdi won't allow me to bring it in when using PythonWin, although it works in the QGIS python command window.

Here is my attempt at trying to import the qgis module:

>>> import qgis
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
ImportError: No module named qgis
>>> from qgis.core
Traceback (  File "<interactive input>", line 1
    from qgis.core
                 ^
SyntaxError: invalid syntax
>>> import qgis.core
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
ImportError: No module named qgis.core

>>> import os
>>> import sys
>>> x = sys.path.append(r'C:/PROGRA~2/QGIS2~1.14/apps/qgis-ltr/./python\\qgis\\__init__.py')
>>> x
>>> os.listdir(x)
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
TypeError: coercing to Unicode: need string or buffer, NoneType foundfrom x

>>> from x import qgis-ltr
Traceback (  File "<interactive input>", line 1
    from x import qgis-ltr
                      ^
SyntaxError: invalid syntax
>>> import x
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
ImportError: No module named x
>>> x = r"C:\Program Files (x86)\QGIS 2.14\apps\qgis-ltr\python\qgis"
>>> sys.path.append(x)
>>> import qgis
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
ImportError: No module named qgis
>>> os.listdir(x)
['analysis', 'core', 'gui', 'networkanalysis', 'PyQt', 'testing', 'user.py', 'utils.py', '_analysis.lib', '_analysis.pyd', '_core.lib', '_core.pyd', '_gui.lib', '_gui.pyd', '_networkanalysis.lib', '_networkanalysis.pyd', '__init__.py']
>>> x = r"C:\Program Files (x86)\QGIS 2.14\apps\qgis-ltr\python\qgis\core"
>>> os.listdir(x)
['contextmanagers.py', '__init__.py']


texts: Create a Polygon on Map using PyQt5. There are no errors but my polygon is not visible on map

texts:
from qgis.gui import *
from qgis.core import *
from PyQt5.QtGui import *
from PyQt5.QtCore import *
import sys
from PyQt5.QtWidgets import QApplication
import os

class MyWnd(QMainWindow):
    def __init__(self):
        QMainWindow.__init__(self)
        QgsApplication.setPrefixPath('D:/SRIJA/Testing/QGis', True)
        QgsApplication.initQgis()
        self.canvas = QgsMapCanvas()
        self.canvas.setCanvasColor(Qt.white)
        elf.lyr = QgsVectorLayer('D:/SRIJA/Testing/QGis/qgis_sample_data/shapefiles/airports.shp', 'airports', 'ogr')
        QgsMapLayerRegistry.instance().addMapLayer(self.lyr)
        self.canvas.setExtent(self.lyr.extent())
        self.canvas.setLayerSet([QgsMapCanvasLayer(self.lyr)])
        self.setCentralWidget(self.canvas)
        actionZoomIn = QAction('Zoom in', self)
        actionZoomOut = QAction('Zoom out', self)
        actionPan = QAction('Pan', self)
        actionPoly = QAction('Polygon', self)
        actionPoly.setCheckable(True)
        self.connect(actionPoly, SIGNAL('triggered()'), self.poly)
        self.toolbar.addAction(actionPoly)
        self.toolPoly = PolyMapTool(self.canvas)
        self.toolPoly.setAction(actionPoly)
        self.poly()
        actionZoomIn.setCheckable(True)
        actionZoomOut.setCheckable(True)
        actionPan.setCheckable(True)
        actionZoomIn.triggered.connect(self.zoomIn)
        actionZoomOut.triggered.connect(self.zoomOut)
        actionPan.triggered.connect(self.pan)
        self.toolbar = self.addToolBar('Canvas actions')
        self.toolbar.addAction(actionZoomIn)
        self.toolbar.addAction(actionZoomOut)
        self.toolbar.addAction(actionPan)
        self.toolPan = QgsMapToolPan(self.canvas)
        self.toolPan.setAction(actionPan)
        self.toolZoomIn = QgsMapToolZoom(self.canvas, False) # false = in
        self.toolZoomIn.setAction(actionZoomIn)
        self.toolZoomOut = QgsMapToolZoom(self.canvas, True) # true = out
        self.toolZoomOut.setAction(actionZoomOut)
        self.pan()

def poly(self):
    self.canvas.setMapTool(self.toolPoly)

def zoomIn(self):
    self.canvas.setMapTool(self.toolZoomIn)

def zoomOut(self):
    self.canvas.setMapTool(self.toolZoomOut)

def pan(self):
    self.canvas.setMapTool(self.toolPan)

class MainApp(QApplication):

def __init__(self):
    QApplication.__init__(self,[],True)
    wdg = MyWnd()
    wdg.show()
    self.exec_()
if __name__ == '__main__':
    import sys
    app = MainApp()

class PolyMapTool(QgsMapToolEmitPoint):

def __init__(self, canvas):
    self.canvas = canvas
    QgsMapToolEmitPoint.__init__(self, self.canvas)
    self.rubberband = QgsRubberBand(self.canvas, QGis.Polygon)
    self.rubberband.setColor(Qt.red)
    self.rubberband.setWidth(1)
    self.point = None
    self.points = []

def canvasPressEvent(self, e):
    self.point = self.toMapCoordinates(e.pos())
    m = QgsVertexMarker(self.canvas)
    m.setCenter(self.point)
    m.setColor(QColor(0,255,0))
    m.setIconSize(5)
    m.setIconType(QgsVertexMarker.ICON_BOX)
    m.setPenWidth(3)
    self.points.append(self.point)
    self.isEmittingPoint = True
    self.showPoly()

def showPoly(self):
    self.rubberband.reset(QGis.Polygon)
    for point in self.points[:-1]:
        self.rubberband.addPoint(point, False)
        self.rubberband.addPoint(self.points[-1], True)
        self.rubberband.show()


texts: Geocoding: Match address is doubling up text

texts:
I'm geocoding some addresses, and using our address data as the locator service's reference. The oddest thing is happening, though, which I've never encountered. The 'Match Address' address is doubling up the street names - giving me a lower Match percentage.
 I'm attaching some photos - if you can give some insight on how I can correct it, it would be a huge help!

As you can see, the Match_Addr is taking the address entered in from Addr field, and doubling it up. It's finding an address and giving me a 100% Score, but keeps the Status as Tied.


texts: Orfeo toolbox and MacOs, installation fails

texts:
I am trying to install orfeo toolbox in my MacAir with Maverics. I am using homebrew but the installation fails with the following error log:

==> Downloading //downloads.sourceforge.net/project/orfeo-toolbox/OTB/OTB-3.20/OTB-3.20
Already downloaded: /Library/Caches/Homebrew/orfeo-3.20.0.tgz
==> cmake .. -DCMAKE_C_FLAGS_RELEASE= -DCMAKE_CXX_FLAGS_RELEASE= -DCMAKE_INSTALL_PREFIX=/usr/

-- Configuring incomplete, errors occurred!
See also "/tmp/orfeo20150720-4244-1cjwrtd/OTB-3.20.0/build/CMakeFiles/CMakeOutput.log".
See also "/tmp/orfeo20150720-4244-1cjwrtd/OTB-3.20.0/build/CMakeFiles/CMakeError.log".

READ THIS: //git.io/brew-troubleshooting

These open issues may also help:
orfeo 5.0.0 (https://github.com/Homebrew/homebrew/pull/41363)
orfeo fails to build on 10.10.3 (https://github.com/Homebrew/homebrew/issues/41275)

I went to the github forum and the issue has been opened, but it was not 100%clear which is the solution


texts: Transforming from PSA56 to NAD83 using ArcGIS Desktop?

texts: I need to transform some data from PSA56 to NAD83.  It doesn't have to be that rigorous as it's purely for a context map.  I'm guessing I will need to go to WGS84 and then to NAD83.  If anyone has any suggestions for best practice please let me know.  Using ArcGIS10


texts: Opening GeoJSON in QGIS: projection problem?

texts:
I'm attempting to use QGIS to convert a GeoJSON to DXF. The problem is that the GeoJSON import is distorted, with the latitudes compressed. This isn't just a problem with the screen's aspect ratio: the distortion carries through to the DXF.

Here's a sample GeoJSON: http://geojson.io/#id=gist:nkoren/8c26a440d9dcff96872b7e7743821803&map=18/51.50123/-0.14192

And this is what it looks like when I import it as a vector layer QGIS 3.6:

My assumption is that this was a projection problem, however, both the layer and the document itself are WGS84, as is the GeoJSON, of course. I'm a QGIS novice, so is there another place that I should be looking?


texts: Slope and hillshade layer not continuous in QGIS

texts: I have been trying to create slope and Hillshade rasters from a DEM raster, which is in UTM11N. I tried different ratios (vertical to horizontal) but everything is in meters so it should work fine. I get these stepping slopes going from 0˚ to 50˚ and back rather than a continuous slope. I was wondering if it was normal to have these lines rather than a continuous slope. If not what can I do to avoid this?


texts: How to use Variable in Select Layer By Attribute in Model Builder

texts:
I have a Variable Types (Multiple Strings) in my model and would like to use the variable with the tool Select Layer By Attribute and a where clause TYPE IN (%Types%).
However, %Types% is a semicolon (;) separated string, e.g. 'Type1;Type2' but what I need is a list of strings (which I would have expected considering that it's a variable with multiple values).
Does ArcGIS Pro has any built-in function to use the variable Types for this use case, or do I have to write a Python script to accomplish this?
I haven't found anything in the SQL reference on the ArcGIS Pro's help page:
https://pro.arcgis.com/en/pro-app/help/mapping/navigation/sql-reference-for-elements-used-in-query-expressions.htm


texts: Looping a location-allocation problem in ArcMap's Model Builder

texts:
I have developed the model below in which I iterate through a series of values for a variable in the Location-Allocation type of problem. The variable is the number of facilities to solve. The model works fine but it is forced to load the facilities and demand points every time it runs.

I cannot figure a work around to use the previously loaded facility candidates and demand points. Since they are always the same this consumes a significant portion of processing time. Does anybody know how to force this variable to change at the "Solve" level?


texts: Using update mode with the SDE30 writer in FME

texts:
I'm recieving the following error when performing an Update operation against a table in ArcSDE using the FME 2011 and the SDE30 writer.

Error retrieving dbtune keyword list:SDE_ERROR_CODE(-66) Invalid parameter value passed to function

Any ideas why this might be?
The table was created using the geodatabase writer in FME 2011.
UPDATE 13/11/2014
I have created a small test harness which is reproducing the error.

If I switch the sde30 writer to:

table exists = yes

then the error goes away,

although the update doesn't work - an insert operation happens instead even with writer in UPDATE mode and fme_db_operation exposed and set to 'UPDATE', however when i switch back to table exists=no, the error comes back.


texts:
I have a route in QGIS that visits every road/edge at least once. By using a marker line I can see the direction of the route, but it doesn't show the order of the segments (see the screenshot below). Is there any style/visualization that will allow me to see the order of the roads in the route?

Currently the route is a single polyline, but I can split it into segments if that will make it easier (the route is generated with Python code).


texts:
I am using r.patch extension in QGIS to join a number of LiDAR tiles, some of which overlap hence i am carrying out this task before utilising in a hydraulic model.

R.Patch works fine, however the precision of the height data is given to 13/14 decimal places - this is causing issues in using the resultant Ascii file in ICM where it cannot read the data.

Is there a solution to reduce the number of decimal points to 2 places?

I am using Version 2.8 as my company hasn't updated the software.


texts:
I have trouble optimizing a query for analytics.

Ok so I have bunch of table (lets say, cable, chamber). All this table have a geometry. This geometry can be linestring, point, etc... (all the same srid). All this table also have a FK to a table project. So we can track project on entities.

Now SQL:

CREATE TABLE cable (
   cab_id uuid primary key default uuid_generate_v4(),
   cab_proj uuid reference project (proj_id),
   the_geom geometry(linestring, 4326)
);

CREATE TABLE chamber (
   chb_id uuid primary key default uuid_generate_v4(),
   chb_proj uuid reference project (proj_id),
   the_geom geometry(point, 4326)
);

CREATE TABLE project (
   proj_id uuid primary key default uuid_generate_v4(),
   proj_desc text
);

What I want is a convex hull around project, dynamically updated, so I have created a view :

CREATE VIEW v_project_limits AS
  SELECT p.proj_id,
    p.proj_desc,
   COALESCE(
     st_collectionextract(
       st_multi(
         st_union(
           st_convexhull(
             st_collect(
               ARRAY [ cab.the_geom, chb.the_geom ]
             )
           )
          )
         ),
       3),
       st_multi(st_geomfromtext('POLYGON((0 0, 0 0, 0 0, 0 0))' :: TEXT, 2154))) AS geom,
  FROM project p
    LEFT JOIN cable cab ON (cab.cab_proj=p.proj_id)
    LEFT JOIN chamber chb ON (chb.chb_proj=p.proj_id)
  GROUP BY p.proj_id, p.proj_desc

On my real data I have 6 tables to analyze. This query makes +4s to execute. I tried query planner explainer like http://tatiyants.com/pev/#/plans but didn't get much result...


texts: How to update an OpenLayers layer from GeoJSON?

texts:
OpenLayers can parse GeoJSON in order to create an array of features. So, it is possible to create a layer from GeoJSON input by simply parsing it and then adding the features to a new vector layer.

So, my current solution sends an AJAX request with the current map extent as a parameter and when the server returns the GeoJSON data, it parses them, then clears the layer and puts the newly received features into the layer.

However, I wonder if OpenLayers can do this automatically. Is it possible to create a layer that requests new GeoJSON data after being moved (e.g. using the BBOX strategy), parses the data and then adds them to the layer? The same way it is done with WFS, for example.

I tried this code:

vectorLayer = new OpenLayers.Layer.Vector("Address points", {
    projection: MAP_PROJECTION,
    strategies: [new OpenLayers.Strategy.BBOX()],
    protocol: new OpenLayers.Protocol.HTTP({
      url: "my-webserver-that-gets-map-extent-and-returns-GeoJSON-encoded-features",
      format: new OpenLayers.Format.GeoJSON(),
      callback: function() { alert('Did at least something happen?'); }
    })
});

Firebug shows that the AJAX request is really being sent every time when I move the map and the server really replies with the GeoJSON data. I used the OpenLayers GeoJSON parser to verify if those data form a correctly structured FeatureCollection, and they do.

However, the features do not appear on the map. I tried the OpenLayers.Control.DrawFeature to draw some points to the map manually and it works ok, so the layer is definitely added and set to visible. The callback function apparently isn't invoked at all (as there is no alert message).

Any hints as to what could be wrong, please?


texts: Two r.mapcalc calculations return a much bigger raster cell size than the input raster

texts:
I tried the following calculations with a GeoTiff sample, having a cell size of 25m:

Calculating the togographic exposure (TOPEX) and
estimating the mean topographic curvature

Both, the resulting raster data from calculation 1. and 2. have a raster cell size of approximately 765m. Even when I choose a 1m GeoTiff, the resulting raster has a cell size of 765m.

Obviously it doesn't depend on the calculation.

Could be the problem that the input file is not type of .dem, .adf or .asc??


texts: Selecting specific featureclasses in result list of arcpy.da.Walk()?

texts:
I am a complete (!) novice in Python.

I have a list of access DB (MDB) that each have 1 or 2 featureclasses that i have to select

example :

B:\map1\DB1.mdb\Feature1
B:\map1\DB1.mdb\Feature1_Corrected
B:\map2\DB2.mdb\Feature1
B:\map2\DB2.mdb\Feature1_Corrected
B:\map3\DB2.mdb\Feature1
B:\map4\DB2.mdb\Feature1
B:\map5\DB5.mdb\Feature1
B:\map5\DB5.mdb\Feature1_Corrected

and so on

unpredictible which Db has 1 or 2 featureclasses

Now I want to append feature1 - if only feature1 exists - or feature1_corrected if both exist to a featureclass in another GDB

I have found the 'walk' function and can append a list to the new GDB but i am a lost in how to write the correct if/then/else statements to select the correct featureclasses during the walk.

I added the script so far below:

# -*- coding: utf-8 -*-

# script to append features that are located in different DB in different maps
# Bart De Wachter - länsstyrelsen Jämtland
# januari 2016
#
# step 1: create list of relevant sjo features - use "walk" -
# step 2: loop for every filename in list
## step 2a: append to feature  (FISK.gdb)
## step 2b: add vattensystem in correct attribute - get info from mapname
# step 3: same procedure for VDR

import datetime

## Timestamp shows
print 'Startar skript vid: \t\t' + str(datetime.datetime.now())
starttid = datetime.datetime.now()

#start with import of relevenat modules
import arcpy, codecs, os
print u'Import av bibliotek färdig: \t' + str(datetime.datetime.now())

# define workspace
workspace = r"B:\GIS\fisk_kunskap\test"
print 'directory = ' + workspace

### STEP 1 ###
# creat list of featureclasses
feature_classes = []

# walk through all directories and databases and find polygon Featureclasses
# works for .MDB and .GDB
walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Polygon")

for dirpath, dirnames, filenames in walk:

        for filename in filenames:
            feature_classes.append(os.path.join(dirpath, filename))

print feature_classes

### STEP 2 ###
# keep only relevant featureclasses
#

### STEP 3 ###
# appends all från feature_classes (list) to specific Featureclass in .GDB
# outfile is where objects are appended
# Schema type = tests if fields are the same; NoTest doesn't check

outfile = r"B:\GIS\fisk_kunskap\Out\OUT.gdb\sjo_all"
schemaType = "NO_TEST"
fieldMappings = ""
subtype = ""

# Process: Append the feature classes into the empty feature class

arcpy.Append_management(feature_classes, outfile, schemaType, fieldMappings, subtype)


texts: Uploading Attachment > 10 MB to ArcGIS Online Hosted Feature Service?

texts:
Has anyone had success uploading an attachment larger than 10 MBs to an ArcGIS Online hosted feature service?

I know that this is not possible via a web browser.  I'm trying to upload a video (.mp4) around 50 MBs.  I broke the video into smaller chunks so that each .mp4 is less than 10 MBs.  I tried the following python script using the uploadPart and addAttachment REST functions to achieve this.  The parts upload successfully, and I then commit and add the attachment to the feature using the commit and addAttachment functions.  However, the attachment is not an .mp4.  It's an application/octet-stream of only 41 bytes.

Here is the python script I'm running:

import requests, json, os

# Disable warnings
requests.packages.urllib3.disable_warnings()

# Variables
username = "agol"
password = "*****"
id = 1252
TempDir = r"D:\Projects\Treddyfrin\Temp"

# Generate AGOL Token
tokenURL = 'https://www.arcgis.com/sharing/rest/generateToken'
params = {'f': 'pjson', 'username': username, 'password': password, 'referer': 'http://www.arcgis.com', 'expiration': 21600}
r = requests.post(tokenURL, data = params, verify=False)
response = json.loads(r.content)
token = response['token']

# Register Item
print("Registering Item")
uploadURL = 'https://services2.arcgis.com/o9a0D2ZwwG8CDFVQ/arcgis/rest/services/Sewers/FeatureServer/uploads/register'
params = {'f': 'pjson', 'token': token, 'itemName': 'attachment1252', 'description': 'Attachment for feature 1252'}
r = requests.post(uploadURL, data = params, verify=False)
response = json.loads(r.content)
itemID = str(response['item']['itemID'])

# Upload Parts of video to Registered Item
def addPart(itemID, token, filename, num):
    '''Add part of the input video using a post request through the request module'''

    print("Adding parts of video")

    addpartURL = 'https://services2.arcgis.com/o9a0D2ZwwG8CDFVQ/arcgis/rest/services/Sewers/FeatureServer/uploads/' + itemID + '/uploadPart'
    params = {'token':token, 'f':'pjson', 'partNumber': str(num)}
    image = {'partFile': (os.path.basename(filename), filename, 'video/mp4')}
    response = requests.post(addpartURL, data = params, files = image)
    print(json.loads(response.text))

# Walk through temp directory to get each video part
for a, b, c in os.walk(TempDir):
    for name in c:
        print(name)
        index = int(name.replace('output', '').replace('.mp4', ''))
        addPart(itemID, token, os.path.join(TempDir,name), index)

# Commit upload
print("Commiting upload")
commitURL = 'https://services2.arcgis.com/o9a0D2ZwwG8CDFVQ/arcgis/rest/services/Sewers/FeatureServer/uploads/' + itemID + '/commit'
params = {'f': 'pjson', 'token': token}
response = requests.post(commitURL, data = params)
print(json.loads(response.text))

# Add attachment to feature
print("Adding attachment to feature")
fsURL = 'https://services2.arcgis.com/o9a0D2ZwwG8CDFVQ/ArcGIS/rest/services/Sewers/FeatureServer/0/' + str(id) + '/addAttachment'
params = {'f': 'pjson', 'token': token, 'uploadId': itemID}
response = requests.post(fsURL, data = params)
print(response.content)


texts: Export Features to Geodatabase Created in same Python Script

texts:
I am new to learning Python for GIS. I am currently trying to create a script which creates a new geodatabase, then using a list, iterates through the features of an existing geodatabase and exports these features to the new geodatabase I just created earlier in the script. I have ended up receiving the error "TEST_GDB (the new created geodatabase) does not exist or is not supported" multiple times. I researched my problem, and got the closest to a solution at this link Using geodatabase feature dataset created in ArcPy script as output/input location later on in same script? however I am still receiving the same error. It is worth noting that existing geodatabase from which I am exporting the data from is a default geodatabase for an ArcGIS Pro project.
import arcpy
# set initial environments
arcpy.env.workspace = r"C:\Users\My_Name\Documents\GIS\ODU_Home\PYTHON_ODU\Python_Script_Tool_Test\Python_Script_Tool_Test.gdb"
arcpy.env.OverwriteOutput = True

# create a list for all feature classes within old or existing geodatabase
OldData = arcpy.ListFeatureClasses()
print("Features Listed")

#Initialize variables to use in CreateFileGDB tool
out_folder_path = r"C:\Users\My_Name\Documents\GIS\ODU_Home\PYTHON_ODU"
out_name = "Test_GDB"

# create a new geodatabase & set as new workspace
GDB = arcpy.CreateFileGDB_management(out_folder_path, out_name,)
print("New .gdb created")
ws = str(str(GDB) + "\\")
arcpy.env.workspace = ws
print("New Workspace")

exportList = []
print("New list initialized")

#export feature classes from old gdb to new gdb
for Old in OldData:
    arcpy.FeatureClassToGeodatabase_conversion(Old, r"C:\Users\My_Name\Documents\GIS\ODU_Home\PYTHON_ODU\Test_GDB")
    exportlist.append(Old + "_NEW")
print("Existing Feature Classes Exported")


texts: Printing error message to console or log file in ArcPy?

texts:
My code is as follows:

Essentially, I have one table in my database that is not able to convert into my FGDB. When I run on command line,  I am able to get an error message of:

  WARNING: Failed to convert Database Connections\DB.sde\Tablexxx. ERROR
  999999: Error executing function.

  Failed to execute (CopyRows).

But i am not able to print/write this error message to console (via print statement) or to a log file (via log.write). It seems that in my code, the except statements are not even being reached, even though an error is occuring?