texts: Remove legend nodes of a categorized layer in a layout legend using PyQgis

texts:
I am trying to find a way to remove specific legend nodes in a layout legend for a vector layer styled with categorized symbology using the PyQgis-3 API.
For example, in a QGIS project I have one categorized layer, a few grouped layers and a number of other ungrouped, vector and raster layers.
How can I remove one or more of the legend nodes belonging to the categorized layer using PyQgis? It is a trivial task to remove them manually via item properties in the layout interface but I can not find any way to accomplish this using PyQgis.
The layout legend (annotated with model rows, groups, parent indices and child node positions) looks like this:

If I access the legend model and get its row count like this:

d = iface.openLayoutDesigners()[0]
l = d.layout()
legend = l.selectedItems()[0]
m = legend.model()
print(m.rowCount())

the row count method returns 8 rows. This tells me that neither the category nodes, nor the child nodes of the group are considered to be separate rows in the model. However the legend nodes of the ungrouped layers and the child nodes of the grouped layers can be removed quite easily.
First an ungrouped layer:

# remove rows with no parent
d = iface.openLayoutDesigners()[0]
l = d.layout()
legend = l.selectedItems()[0]
m = legend.model()
legend.setAutoUpdateModel(False)
m.removeRows(4, 1) # pass row index & number of rows to remove
legend.refresh()

Then a member of the ‘Collecting Targets’ group:

# remove rows from a group
d = iface.openLayoutDesigners()[0]
l = d.layout()
legend = l.selectedItems()[0]
m = legend.model()
idx = m.index(1, 0) # get index of parent
legend.setAutoUpdateModel(False)
m.removeRows(1, 1, idx) # pass position, number of rows to remove & index of parent
legend.refresh()

However, if I use the code above but get the index of the categorized layer and pass that to the removeRows() method as the 3rd argument (parent), no rows are removed. There is no error encountered- it just does not work.

I can use the method legendRootRowCount() to return each category of the categorized layer as a row.

I can also access the nodes of the categorized layer and store the parent and its nodes in a dictionary like so:

d = iface.openLayoutDesigners()[0]
l = d.layout()
legend = l.selectedItems()[0]
layer_and_nodes = dict()
m = legend.model()
idx = m.index(0, 0)
n = m.index2node(idx)
nodes = m.layerLegendNodes(n)
layer_and_nodes[n.name()] = nodes
print(layer_and_nodes)

This returns:

{'fs2019shp': [<qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809CFE58>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809CFEE8>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809CFF78>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809E2048>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809E20D8>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809E2168>, <qgis._core.QgsLayerTreeModelLegendNode object at 0x00000013809E21F8>]}

But I can’t work out any way of removing or filtering those QgsLayerTreeModelLegendNode objects.

I tried:

llns = [v for k, v in layer_and_nodes.items()][0]
m.filterLegendNodes(llns)

I suspect I am not implementing this method correctly as this also fails with no errors.

I try to force a refresh of layer nodes for the categorized layer after both the removeRows() and the filterLegendNodes() calls using m.refreshLayerLegend(n) which does not have any effect.

I have scoured the API documentation and searched GIS.StackExchange for similar questions. My research also introduced me to the concept of QT model/ view programming (of which I admit I was completely ignorant) but I have not been able to find a solution.

How can I remove a legend node of a categorized layer from a QgsLegendModel using PyQgis?


texts: Using output feature class field values in SQL query in model builder

texts:
I have two processes in model builder:

Iterate over a "ID_SPATIAL" field values = creating a subset of original feature class (and agregating using Frequency tool)
Creating a different subset using a "Code" field from the original feature class by SQL query

What I need is to use the values from "Code" field within the first subset as an input to SQL query. Model looks like this:

My SQL query could look like this:

SELECT * FROM INPUT2 WHERE Code IN (SELECT Code FROM %ID_freq list%)

The attribute table of objects looks like this:

I have spent last day wondering about this problem. How can I transfer values from field of another process into SQL query?

I don't have a solution yet... but I succeded to make a procedure using a submodels with Merge tool and Collect Values tool of model builder. This is very programming nasty, but since I haven't found a solution it helped.
THe problem is, that this is very very time consuming procedure, and I have to run it over and over again after an hour when ArcMap depletes its inner memory.

Submodel: for each ID_SPATIAL this iterates "Code" field and back query the data (A02_108_not_null_codes) "Code" field for each one. Afterwards, it passes selected features to Merge tool (overall model).

Overall model: merge the results in one layer

This should serve as a typo finder. Once in a map, I can see if some codes were mistyped because these features lay apart from the rest of particular ID_SPATIAL region:


texts:
I am checking effects of reprojecting planar coordinates from one CRS to another and back, keeping the same DATUM, and I was trying to determine the cause of the loss of precision.

I know that deformation grows to infinite values, and that precision will be lost eventually, but I was surprised to see differences in nearby zones. Am I correct to state that it is due to coordinates projected to far-out zones to grow in number size thus limiting the number of decimal places available in double precision arithmetics?

My version of Postgis and linked libraries: POSTGIS="2.3.3 r15473" GEOS="3.4.2-CAPI-1.8.2 r3921" SFCGAL="1.3.0" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER

Below a reproducible example:

select *   from (SELECT 'original f32' as crs, ST_x(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)')) as forward, ST_x(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)')) as back
UNION
SELECT 'original f32->f33->f32' as crs, ST_x(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32633)) as forward,
ST_x(st_transform(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32633),32632)) as back
UNION
SELECT 'original f32->f34->f32' as crs, ST_x(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32634)) as forward,
ST_x(st_transform(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32634),32632)) as back
UNION
SELECT 'original f32->f35->f32' as crs, ST_x(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32635)) as forward,
ST_x(st_transform(st_transform(
ST_GeomFromEWKT('SRID=32632;POINT(500000.0 4000000.00000 0)'),32635),32632)) as back
) a order by crs


texts: How to create a missing dbf file for shapefile?

texts:
I have a shape file (shp) which is missing its dbf. I tried to load this into Postgis using shp2pgsql and got the error.

  could not open shapename.dbf

I opened it in QGIS, without issue, and when I clicked Open attribute table, I could see a column based on the shp's gid, but could not use the field editor to add a new field. Instead, I used, Save as and saved as a new shp file.

This new shp file, once added to the map, did allow me to add fields and I can now successfully load it into Postgis using shp2pgsql.

Question: while the above works, it is a bit of a faff, and I am wondering if I have missed some obvious trick, either in shp2pgsql, QGis, ogr2ogr or some other application not mentioned.


texts: Removing Overlapping Outlines Within a Single Feature Class

texts:
I have a polygon feature that is overlapping itself. Both of the polygons in the image below are from that same feature. I do not want to dissolve the boundaries, but I only want one of the outlines to present. They also have different line styles for their outlines, so I want to be able to choose which one presents depending on the contact.

Ex. - Polgyon A overlaps Polygon B and Polygon C. At the contacts, I want Polygon A's outline to present at contact with Polygon B, but at Polygon A's contact with Polygon C, I want Polygon C's outline to present.


texts: Edit UI - Python Init function localization with QTranslator?

texts:
I have  interaction with layer defined by custom designed UIs.
In layer properties the  following settings: Edit UI - customline.ui and Python Init function: customline_set.py.

#coding: utf8
from PyQt4.QtCore import *
from PyQt4.QtGui import *

nameField = None
myDialog = None

def formOpen(dialog,layerid,featureid):
    global myDialog
    myDialog = dialog
    global nameField

    translator = QTranslator()
    loc = QSettings().value('locale/userLocale')[0:2]
    print loc
    filename = "\\i18n\\customline_%s.qm" % loc
    print filename
    translator.load(filename)
    QCoreApplication.installTranslator(translator)

    print QCoreApplication.translate('editForm',"Generals")
    # return translated string ok

The code work without error but the UI it's not displayed in selected language.

Can anyone point me in the direction of the mistake here?


texts: Exporting data from ArcGIS Desktop to .txt?

texts:
I have a file-geodatabase.  It is a road network of a city.

Now I want to export the data into a .txt file for map matching.

I have installed ArcGIS Desktop and not used that software before.

What should I do next?

open the data with ArcMap

txt format

I am completely new to GIS field!


texts:
I am new with QGIS.

I want to study the performance of different GPS receivers, and to do so I need to compute the error made by each one. The problem is that I don't have any reference to be compared with the coordinates provided by the receivers (which are installed on a car in movement), so I want to generate it myself using a GIS software.

The idea is to stick my GPS waypoints to the closest road (just as Google Maps app) to create a new route that would remain always within the road. And then I would like to (somehow) import this new data to matlab (for example, in CSV format) in order to compute the positioning error.

Is it possible or I should try something else?

I have seen that there are many functions available using PostGIS, and I have tried with ST_ClosestPoint (http://postgis.net/docs/ST_ClosestPoint.html) without any luck.

What is the correct way of using this function if I want to use the polylines of my city and the position fixes of my receiver as inputs ?


texts: Python Script to Add Fields to Feature Classes

texts: I have to try to add a field called 'Name' to a bunch of different feature classes using python, and then access the name of the feature class to fill in the field. I have no idea how to even start this! I know I need a loop but that's all I know.


texts:
I am looking for a way to calculate the segment lengths of a line network. I am probably missing just the right vocabulary to search for the best approach, so I guess asking here would be a good idea. I have prepared an image of what I want to achive:

Please note that I am not interested in the length of each line, I know how to calculate that. What I want to calculate is the length between the squares, the intersections that is. Bonus would be to have the lengths towards the end points of a line.

I imagine that pgRouting is able to do that, but I'd prefer a solution without Postgis.

Is this already built into QGIS with the Network analysis library?

I am unsure what to look for.


texts:
I'd like to only display an inset map, inset map scale bar, frame and leader if an associated layer is displayed. Similar to Associate Point Layer, Annotation Group and Line Layer in Map Legend but for an inset map rather than annotation.

Any thoughts how to approach this?


texts: Determine values of a point based nearest point values - QGIS

texts:
I have some weather stations that receive and publish precipitation data every day.
Then I have polygons (can be center points if necessary) that correspond to interest areas.

I want to estimate the precipitation for each area with the stations within a maximum buffer of the center point.

It's not possible to use the default interpolation tools because the stations not always surround completely the polygons.

Any idea to solve this problem?


texts: Can Grid Index Features in ArcGIS Desktop fit over oddly shaped polygons?

texts:
I am using ArcGIS 10.0 where I work.

One of the clients is wanting a map book as well as a locator map that shows where each page of the map is showing.  In order to do this I must first create an Index Grid with the Grid Index Feature tool.  Where I am running into trouble is the features associated with the shapefile are oddly shaped and vary in size.

Is there a way to make the grids fit over each feature in the shapefile?

Currently the Grid Index Feature tool only wants to give me boxes that fit over multiple and partial shapes.


texts: OSM in QuickMapServices not displaying baselayer image on QGIS

texts: [I am attempting to load a QuickMapServices OSM baselayer image on QGIS 3.4 LTR and nothing is displayed, even though it says that the layer has been added. Nothing seems wrong with my plug-ins and I cannot seem to figure out why.


texts: How to display changes in lake's surface area over time in a polygon bounded region, in a table in earth engine?

texts:
The problem statement is that a region of interest is given.

I need to find all the lakes in a polygon bounded region using the NDWI index for water bodies, which are at a height of more than 1500m. Then display the changes in the area of lake's surface water starting from the year 1984 till 2018 on a 2-year interval in a table like structure in Google Earth Engine. I have used Landsat 5 and 7 data.

I need to display the results in the polygon marked region in a table sort of structure in the following format:- Rows - (Lake 1, Lake 2, Lake 3... Lake n) Columns - (Surface Area in 1984, Surface Area in 1986, ....2018) - surface area in square km.

How should I go about doing it?

I have created the following code:

    var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR"),
    roi = /* color: #d63000 */ee.Geometry.Polygon(
        [[[77.68018052033392, 31.0091863423649],
          [80.36084458283392, 29.5863837026937],
          [80.98706528595892, 30.233944725529557],
          [78.85571762970892, 31.57247908795104]]]),
    elevation = ee.Image("USGS/SRTMGL1_003"),
    l7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_RT_TOA");

// Defining a filter that gives results for months September to November
var mf = ee.Filter.calendarRange(9, 11, 'month');

Map.centerObject(roi, 9);

// Definig the image collections
var img1 = ee.ImageCollection(l5
            .filterDate('1984-01-01','1999-12-31')
            .filterBounds(roi)
            .filter(mf));
var img2 = ee.ImageCollection(l7
            .filterDate('2000-01-01','2018-12-31')
            .filterBounds(roi)
            .filter(mf));

// add NDWI band function
var addNDWI = function(image) {
  return image
    // NDWI
    .addBands(image.normalizedDifference(['B2', 'B4']).rename('NDWI'))
};

// map add ndwi function to both image collections
var image1 = img1.map(addNDWI);
var image2 = img2.map(addNDWI);

var test2 = image2.median();
Map.addLayer(test2.clip(roi), {}, 'image2');

// masking function for elevation
var mask = function(image) {
  var elevMask = elevation.updateMask(elevation.gt(1500));
  return image.updateMask(elevMask);
};

// masking function for NDWI
var mask1 = function(image){
  var ndwiMask = image.select('NDWI').gte(0.4);
  return image.updateMask(ndwiMask);
};

// map masking functions
var maskedImg = image1.map(mask).map(mask1);     //For 1984-1999
var maskedImg2 = image2.map(mask).map(mask1);    //For 2000-2018

// getting rid of bad pixels
function trimImage(image){
  return image.clip(image.geometry().buffer(-5000));
}

maskedImg2 = maskedImg2.map(trimImage);
print(maskedImg);
print(maskedImg2);

//set up params to view NDWI band
var ndwiViz = {bands: 'NDWI',min: 0.3, palette: ['00FFFF', '0000FF']};

var test = maskedImg.median();
var test1 = maskedImg2.median();

Map.addLayer(test.clip(roi), ndwiViz, 'NDWI 1984-1999');
Map.addLayer(test1.clip(roi), ndwiViz, 'NDWI 2000-2018');

Also I'm trying to get rid of the white patches in the output.

This is the Earth Engine link to the same


texts: Appending the features to PostGIS with PyQGIS

texts:
I have two layers with same schema and same name but its contain different features.
I need to append those features with PostgreSQL db layer. All three Layers(db layer) schema and name is same .
def Transfer_data_to_Postgis(Tablayer):
  uri = QgsDataSourceUri()
  uri.setConnection(self.Hostname, self.Portid,self.Dbname,self.Username, self.Password)
  uri.setDataSource(self.Schema,Tablayer.name(),"geom","Version=1")
  dblayer= QgsVectorLayer(uri.uri(),Tablayer.name(), "postgres")
  dblayer.startEditing()
  print (dblayer.isValid())
  print (dblayer.isEditable())
  for fea in Tablayer.getFeatures():
      dblayer.addFeature(fea)
  dblayer.commitChanges()

Layers=QgsProject.instance().mapLayers().values()
for tablay in Layers:
   Transfer_data_to_Postgis(tablay)

I tried above code its pushing single layer for other layer script working fine but features not appending.
What's wrong with this script and how to append features  to PostGIS db layer?


texts: ERROR 000210 Cannot create output 'name' Failed to execute CopyFeatures

texts:
I'm learning ArcPy, when I run this code from Python on ArcMap it works. But it does not work from the separate Python shell.

import arcpy
from arcpy import env
import os
output_dir = "C:/Data"
env.workspace = os.path.join(output_dir, "Output.gdb")
env.overwriteOutput = True

output_features = "name"
origins = "C:\Users\username\AppData\Roaming\ESRI\Desktop10.5\ArcCatalog\Connection to sde.sde\SDE.tablename"

arcpy.management.CopyFeatures(origins, output_features)

I'm getting arcgisscripting.ExecuteError: ERROR 000210: Cannot create output 'name' Failed to execute (CopyFeatures).


texts: Providing Enterprise License Server Management?

texts:
I was hoping to get some information from others who have deployed an enterprise GIS License Manager for ArcGIS 10.x software.  Honestly, I would like to see if our current setup is adequate or if we should venture into a new setup architecture.

Currently, we have multiple departments who each purchase their own licenses.  To ensure departments only have access to their purchased licenses, ArcGIS License Managers have been installed for each department and their licenses authorized on each LM respectively.  Some of these LM's are on their own dedicated virtual server, and some are installed on client machines in the department (not ideal and we plan to change this immediately).  Our idea to both keep the licenses separated, and also not use any of our farm resources is to place physical boxes (ones that have been replaced via PC upgrades, etc.) in each department, install LM on the box and authorize the licenses.  Each department would have it's own box, and we can administer these via remove access.

Here are questions I have:

Is it at all possible to run multiple instances of LM on one machine, and just change the outgoing port so each instance has it's own dedicated I/O access? (i.e. Dept A uses 27000, Dept B uses 27001, etc.)
I have heard (and read somewhere), that individual licenses on LM's can be managed down to the user level via changing the LM configuration file.  If this is so, have any of you implemented this setup and how successful has it proved to be?

Again, we are just wanting to ensure we are implementing a practical architecture.  We've used the above setup for years so I believe it's time we revisit this and see if better practices exist.

Here are some of the posts/articles we have looked over for ideas so far:

Managing ArcGIS floating licenses
hardcoding port number is ArcGIS License Manager


texts:
I have 64 bit as lat/long (32 bit lat , 32 bit long). I don't know its original format (and type) but I know its position on Google map. I extract lat long as below (in MATLAB):

lat=uint32*180/2^31
long=uint32*180/2^31

but when I see lat long on the Google map I have many degree error in lat long. How can I convert it with no error? Does any one know this lat long format (or type)?


texts:
I have a feature class. I want to edit the alias. When I go into the Properties of the FC and try to edit the alias, nothing happens. The field will not accept focus and changes are not allowed. How do I change the alias of a feature class? This could be done via code, ArcSDE command line, or through the ArcMap / ArcCatalog UIs.

This is a feature class in an enterprise ArcSDE (Oracle).
The user is the schema owner.
Attempted to change the alias using ArcMap and ArcCatalog.


texts: Linked layers / virtually merge layers in QGIS?

texts:
I have 20 different regions (let’s say 20 provinces). For each province I have let’s say 50 different attribute tables / layers containing data for the province (roads, bridges, towers, …). I would like to show the Data for all provinces in one QGIS file. So I have to import all these layers into one QGIS file ending up with 20 road layers, 20 bridge layers, 20 tower layers,….
This brings up 2 problems:

Due to the amount of layers it is extremely confusing to work with the file even if I group them.

When I change the symbology of the tower layer I have to change it 20 times, right?

So what I did was I merged the road data, bridge data, tower data into one attribute table each, ending up with one layer containing the roads of all provinces, the bridges of all provinces, the towers of all provinces. Working with the file now is more easy due to the extremely reduced number of layers and when changing layer symbology e.g. for towers I only have to do this once.
The problem is though that when I get updated Information on a province I have to manually merge the data with each attribute table again (new infos on roads, new info on towers ….), deleting the former information in the related tables which all together takes a lot of time and is confusing and error-prone.
I hope this was so far understandable.
What could be a solution is to put all the data (shape-files) in a folder on the harddrive and just virtually merge them in QGIS, like they were just linked. They would show up as one layer in QGIS but the referred data lies separated on the drive. When there is new data available for one province I would then just replace the files in the folder on the harddrive and next time QGIS starts, the virtually merged / combined / linked table would contain the updated information.
Is there a way to do this in QGIS?
Or is there another way of handling this? One way would be to have separate QGIS files for each province which would be my dream anyway, since working with the file is really slow, and saving the file takes ages. But the problem with this solution is that we quite frequently have to change layer symbology. Each time we change the line weight for roads you would have to open all files one after the other and change the line weight of the road layer to be consistent (at least as far as I know there is no other way). This is not acceptable. Does anyone have a solution following this path?


texts: I created text diagrams for my polygons and want to juse them as labels in the print layout. Is there  way how to adjust the Textsize of the diagrams? I was thinking about labeling the centriod of the polygon with .svg, but how can export the diagrams as SVG?


texts: My doubt is that I have to do a calculation between 4 rasters. Each raster is classified in 1 to 5. So the final calculation must be a numerical sequence, like: 1111, 1321, 1555, 1432.... for example. Anyone knows how can I do that? Can be in ArcGis or QGIS?


texts: rasterio driver - weird bug that throws no error

texts:
I am using rasterio to convert some png files to geotiff files.
My input is a geotiff and I do some processing to output a geotiff file using the following command:
        temp = cv2.imread(file)
        temp = cv2.cvtColor(temp, cv2.COLOR_BGR2RGB)
        src = rasterio.open(geo_file)
        tmp_outpath = os.path.join(file.replace('.png', '.tif'))
        with rasterio.open(tmp_outpath, 'w', driver='GTiff', height=src.height, width=src.width,
                           count=temp.shape[2], crs=src.crs, transform=src.transform, dtype=rasterio.uint8) as dest_tiff:
            for i in range(temp.shape[2]):
                dest_tiff.write(temp[:, :, i], i + 1)

where file is my processed png file and has the same dimension as the geo_file.
In the remote machine, the geotiff is not saved with the georeferencing I had extracted out from geo_file. The bug is occuring in a remote VM but works perfectly fine on my machine.


texts:
I'm using OpenLayers 4.6.5. I have a VectorLayer with a couple of circle features that have the same zIndex. The problem is, that the features texts are overlapping although I've set the the opacity of the text and the fill to 1.0. The following image shows two circle features with texts '3' and '4' overlapping. I want that one circle is entirely below the other. It doesn't matter which one but I cannot change the zIndex.

var markerStyle = (feature: Feature) => new Style({
    text: new Text({
      text: feature.get('text'),
      font: 'Normal 8px Verdana, Geneva, Tahoma, sans-serif',
      fill: new Fill({color: 'rgba(255, 255, 255, 1.0)'})
    }),
    image: new Circle({
      radius: 10,
      fill: new Fill({color: rgba(193, 0, 31, 1.0)}),
      stroke: new Stroke({color: 'rgba(255, 255, 255, 1.0)', width: 1}),
    }),
    zIndex: 10
  });

So basically the two features have the same style but different texts. Interestingly enough I have other circle features that overlap only partially and the stroke color of one circle is always below the other circle (i.e. not visible), but the text is not. I cannot change the zIndex as there are too many different features.


texts: passing geocode results object to another file

texts:
I'm using a geocoding plugin (Leafet.GeoSearch) on a map and I'd like the results, being the LatLng point coordinates, etc available as a parameter for GeoServer cql filtered calls.

It doesn't seem to be as easy as I thought it would be. I've read several forums which claim data can't be persistent between pages/files. Could someone enlighten me and or make a suggestion please. I'm in over my head.

here is a snippet of the code that has the array I'd like available

ParseJSON : function(data) {
        if (data.resourceSets.length == 0 || data.resourceSets[0].resources.length == 0)
            return [];

        var results = [];
        for (var i = 0; i < data.resourceSets[0].resources.length; i++)
            results.push(new L.GeoSearch.Result(
                data.resourceSets[0].resources[i].point.coordinates[1],
                data.resourceSets[0].resources[i].point.coordinates[0],
                data.resourceSets[0].resources[i].address.formattedAddress));
        return results;
    }


texts: OGR2OGR import to SQL from geometry to geography

texts:
I'm trying to import a shapefile to SQL Server 2019. This command works:
ogr2ogr -f "MSSQLSpatial" "MSSQL:server=[myserver];database=test;tables=CAPZONE;UID=[myuser];PWD=[mypassword];driver={SQL Server}" "[myshapefile]" -lco GEOMETRY_NAME=GEOM -lco GEOM_TYPE=GEOMETRY -nln "CAPZONE" -a_srs "EPSG:4326" -overwrite -progress -skipfailures -lco UPLOAD_GEOM_FORMAT=wkt

But it creates the column with the data in SQL of type GEOMETRY if I change the parameter GEOM_TYPE=GEOMETRY to GEOM_TYPE=GEOGRAPHY it also works (as in I don't get any errors) but the difference is that with the GEOMETRY type the shapes are correct when I check in SQL instead with the GEOGRAPHY type I only get a full page of color in the "spatial results" pane.
I need the GEOGRAPHY type to be able to run some SQL functions not available with the GEOMETRY type. I've also tried importing the GEOMETRY type and then using some SQL functions to convert to GEOGRAPHYbut the result is the same.
this is my [spatial_ref_sys] table:

Any ideas?
-----UPDATE------------
Adding parameter
-nlt "LINESTRING" or "MULTILINESTRING" or "COMPOUNDCURVE" finally draws the right shapes but they aren't filled with colour for some reason:

anyone know why it isn't filling the up?


texts: How to join number of point features into attribute of nearest line segment using ArcGIS Desktop?

texts:
I have point shapefile with lot of attributes, that represent road accidents and line shapefile with highways. What I want is to join attribute into line shapefile with value, that represent number of accidents which happened on each road segment (line segment in shapefile).  First idea was to create buffer around line segments which was used to locate accidents, but I don't know what to do next.

I need this to create simple line diagram map (exactly this one http://athena.zcu.cz/kurzy/gish/000/HTML/37/Ob98.jpg) and I want to know some possibilities of  cartographic visualization from these data, not only joining number of accidents, but some idea to use attributes in data too. I'm talking about simple line diagram map.

I am using ArcGIS Desktop 10.


texts: I cannot find the script gdalcompare.py within the osGeo distribution of GDAL by http://www.osgeo.org/. It is the 1.10.1 though. Is there any way of downloading that script separately from any other source?


texts:
I have the following activity called activity_esri_map:

<LinearLayout xmlns:android="http://schemas.android.com/apk/res/android"
xmlns:tools="http://schemas.android.com/tools"
android:layout_width="match_parent"
android:layout_height="match_parent"
tools:context="com.example.stevecortis.wmsdemo.MainActivity">

<com.esri.android.map.MapView
    android:id="@+id/arcMap"
    android:layout_width="fill_parent"
    android:layout_height="fill_parent"
    mapoptions.MapType="Topo"
    mapoptions.ZoomLevel="13"
    mapoptions.center="33.666354, -117.903557" />

This is called from the MainActivity:

final MapOptions mStreetsBasemap = new MapOptions(MapOptions.MapType.STREETS);

private MapView arcMapView;

@Override
public void onCreate(Bundle savedInstanceState) {
    super.onCreate(savedInstanceState);
    setContentView(R.layout.activity_esri_map);
    arcMapView = (MapView) findViewById(R.id.arcMap);
    arcMapView.setMapOptions(mStreetsBasemap);
}

From the MainActivity I haven't been able to change anything from tha map. It seems like the findViewById(R.id.arcMap) It is not being assigned to the arcMapView.

I have tried to change other things such as Scale, and nothing seems to affect the actual MapView.


texts: Interpolating points observed with different time duration

texts:
I am doing spatial interpolation with a group of points.

Some points have complete observation of events:

for a period of 24 hours we have observed N events in each point. These points have N events.

While others are incomplete:

for a period of 6 hours we have observed M events in each point. This means that these points have >= M events, since we don't have 24h of observation we just know what is the value for the minimum amount of events.

Although the second group of points have incomplete observations, they are still important for the interpolation. I would like to use them in the spatial interpolation (kriging or co-kriging) with the condition that these are not complete observations.

I would like to interpolate new points and use both complete (24h) and incomplete (6h) points. How can I take into consideration, in a spatial interpolation, points that are incomplete?

How can this be done?


texts: Insert shapefile features into PostgreSQL database using java EE (pragmatically)

texts:
I am working with PostgreSQL 9.4 and Java EE. I want to insert shapefile or layer's features into PostgreSQL, is there any solution so I can insert the whole shapefile into the database? OS is Windows 10.

I have successfully inserted one point into my database (PostgreSQL) table but I want all the features like: polygon, multi-polygon, marker etc.

Whatever I draw on the map should be stored into particular table of DB. I am getting all the features but unable to store them. How I can do this?

    try
    {   Class.forName("org.postgresql.Driver");
        conn=DriverManager.getConnection(url,username, password);
        sql=conn.createStatement();
        result1 = "insert into test(feature) values(ST_GeomFromText("+data+"))";
        result1 = "insert into test(feature,lat,lan) values(ST_GeomFromText('"+feature+"),"+lat+","+lan+"')";
        result1 = "insert into test(feature) values(ST_GeomFromText('POINT("+lat+" "+lan+")', 4326))";
        result1 = "insert into test(feature) values(ST_GeomFromText('"+feature+"("+lat+" "+lan+")'))";
        sql.executeUpdate(result1);
        out.println("data inserted...");

        String nextJSP = "/index.jsp";
        RequestDispatcher dispatcher = getServletContext().getRequestDispatcher(nextJSP);
        dispatcher.forward(request,response);

        out.println("forwarding...");
        //System.out.println(fst);
    }
    catch(SQLException | ClassNotFoundException e)
    {
        e.printStackTrace();
    }
}

Getting table data and project it on map:

<%
JSONObject GJson = (JSONObject) request.getAttribute("geoJson");
//String GJson = (String)request.getAttribute("geoJson");
System.out.println("data at jsp : "+GJson);
%>


texts: Proper syntax for Select Feature Layer by Location

texts:
Can someone let me know how I can use a single feature from one layer as the "Select Features" when using the Make Feature Layer by Location tool with another layer as the input?

In the code below I have no idea what to put in place of BLAH to get this to run.  Suggestions?

import arcpy

Practice_Sales1 = "Practice_Sales1"
Practice_Stops = "Practice_Stops"

import arcpy
fc = "y:/Practice_Sales1.shp"
field = "BOROUGH"
cursor = arcpy.SearchCursor(fc)

# Iterate through the rows in the cursor
for row in cursor:

    #Create Layer from the current row in the Practice_Sales1 file
     arcpy.MakeFeatureLayer_management("Practice_Sales.shp", "Sales_lyr", BLAH)

    # Select all features in Practice_Stops within 1000 feet of the current row in Practice_Sales1
    arcpy.SelectLayerByLocation_management(Practice_Stops, "WITHIN_A_DISTANCE", "Sales_lyr", "1000 Feet", "NEW_SELECTION", "NOT_INVERT")

    # Returns the # of selected features
    arcpy.GetCount_management(Practice_Stops)

    # Replace row value in COUNT field with
    row.setValue("COUNT", arcpy.GetCount_management(Practice_Stops))


texts:
Using QGIS 2.18.16 I am looking at "renaming" or changing the attributes of points to a unique attribute based on the direction of the points.

For example A1, A2, A3, then next row starts over as B1, B2, B3...

I would like to start at the furthermost point SW and gong in a straight line NE to the end until there is no longer points, then start again with the next row to the west.

Any help would be great.


texts: Is there maximum number of .DBFs that can be joined using ModelBuilder?

texts: Is there a maximum number of .DBFs that can be joined using the Join Tool in ModelBuilder? I tried joining 480 .DBFs but results generated only 250 fields in the input table. Is this a glitch? I am using ArcGIS for Desktop 10.2.


texts: Using ArcGIS Pro raster calculator to give cell value based on cells in its neighborhood?

texts:
I have a land cover dataset, I want to find forested cells that are adjacent (touch) agricultural cells.

Is there a way I can achieve this with the ArcGIS Pro raster calculator?


texts: ArcToolbox does not open from ArcCatalog?

texts:
From ArcCatalog 10.3.1 I can´t open ArcToolbox.

Software just freezes and crashes when I press anything.

Does anyone know how to solve this?


texts: Fixing "The proj4 definition of Generated CRS ... is not valid" in custom CRS dialog

texts:
I just installed QGIS v2.18.2 for Mac and can no longer define a custom CRS. I get the following error. I've tried using the same parameters as existing CRS and still get the error.

Ive tried complete reinstall and restart but am still having this problem. Is there anything else I could try?


texts: Seeking shapefiles for Marshlands, swampy areas in UK?

texts:
Can anyone point me in the direction of shapefile data for all the coastal and inland areas of saltmarsh, wetland and inland marshes for the United Kingdom, including Isle of Man?

OS Open Data does not show any wetland/marshy ground on any of it's products!


texts: How to force fixed position of lines and gaps pattern in SLD/YLSD with stroke-dasharray or stroke-graphic

texts:
Any kind of dashed effects (simple line or svg/image symbols on line) are rendered correctly at fixed position if the whole line or its start point are currently visible. If we move the map view that the start point of line is outside of it, the line pattern will be shifted like this:

name: LINE
title: LINE
feature-styles:
- rules:
  - name: LINE
    symbolizers:
    - line:
        uom: metre
        stroke-color: '#000000'
        stroke-width: 0.5
        stroke-linejoin: bevel
        stroke-linecap: butt
        stroke-dasharray: 0.75 0.25

On first picture the line start point is bottom left and the line is all visible on map. Its starts from 0.75m dash and after that 0.25m gap and so on. At the of the line end top-right it ends with 0.25 gap. If I move the map a little bit down so the start of my line isn't visible, the end of my line is not ending with a 0.25 gap (picture 2, 3).

Is there any way to force fixed position of stroke-dasharray or stroke-graphic pattern?

Below SVG example: 1 whole line visible (correct); 2,3 part of line visible (incorrect/shifted).


texts: How to combine historic OS sheets in order to create density maps of historic streets

texts:
I really need some advice on creating historical street maps using GIS software (I have access to QGIS and Diva).

Basically, I want to upload a number of historic OS sheets, which I have downloaded from Digimaps in tiff format.

I then want to create dot density (or similar) maps focusing on specific districts, which will illustrate instances of prostitution. Many of the streets I want to map do not exist any more, hence the need for historic OS maps.

I've been fiddling around with the two software packages but I am a complete novice and cannot find out how to do what I need to do.

Could anyone help?

Many thanks in advance!

(Note that this question is not a duplicate as I am looking to stitch together historic OS maps - see my comments in the discussion below).


texts: Raster alignment constraint violated raster2pgsql

texts:
I have a problem when loading multiple raster files into a Postgis database.
When I run this command in the cl

raster2pgsql -I -C -e -Y -F -s 4326 -t 100x100 *.tif public.newTable | psql -d db

the following warnings are returned:

NOTICE:  The rasters (pixel corner coordinates) are not aligned
NOTICE:  Unable to add constraint: enforce_same_alignment_rast

The rasters seem not to have exactly the same alignment. When looking at the corner coordinates, I see that they differ. Two corners that should collapse have the following coordinates (WGS):

Lower Left  (   0.0000000,   9.9996668) (  0d 0' 0.01"E,  9d59'58.80"N)
Upper Left  (   0.0000000,  10.0000001) (  0d 0' 0.01"E, 10d 0' 0.00"N)

I also tried to load the files to R using the raster package, but it didn't work either. However when loading the rasters directly into QGIS there is no problem.

Does anybody know how to resolve this issue? And why does it work in QGIS, but not in R and PostGIS?


texts: How to print automatically Map themes in Label's Map Composer?

texts:
This expression [%@project_path%] prints the path of the qgis project.

How the Map Theme used in the layout can also be printed automatically ?

According to QGIS documentation, a Map Theme is a snapshot of the current map legend that records.

How to reproduce the issue

The Map Theme is set according to these steps :

View tab > Panels > Layers

Click on the Manage Map Themes icon > Add Theme

Give the new Map Theme a name

Open a new Print Layout > click on Adds a new Label to the Layout > create one

Issue

Does a variable or an expression exist to automatically print the name of the Map Theme in the Print Composer?

Answer 1 (not entirely satisfactory)

create manually a Map Theme
open Python console > iface.mapCanvas().setTheme('your theme')
open Print Composer > add a new map to the layout > Item propriety tab > Item ID > write what you want
add a new Label to the Layout > item propriety > click insert an expression button > function editor tab > click '+' button > insert the code below :

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args='auto', group='Custom')
def print_theme(value1, feature, parent):
    return iface.mapCanvas().theme()
Expression tab > print_theme('_your_map_id_')


texts: Thinning LiDAR points to remove redundant neighbors?

texts:
I'm trying to develop a workflow in ArcGIS with 3D Analyst or Spatial Analyst to intelligently thin a grid of spot elevation points created from a LiDAR DEM. The current point spacing is 1 meter. Some CAD software packages have tools to do this kind of thinning directly on a TIN surface, but I would like to port that functionality into ArcGIS. I'm guessing the workflow would involve some kind of nearest neighbor analysis.

What I would like to do is intelligently thin a grid of points to remove points in areas of redundant elevations. For example, a large parking lot may contain 200 points whose elevations are all within 0.25' of each other. This parking lot's elevations could instead be accurately represented in the surface with maybe 10 points.

The workflow would aim to eliminate neighboring points with similar elevations to reduce the overall number of points in the spot elevation grid. User controls would decide what the height and distance thresholds are to determine whether to remove points. Following the thinning, these spot elevations would get fed into a TIN.


texts: Randf in correlation to other field QGIS

texts:
I'm looking for solution that will allow me to create random values with correlation to other field.
For example:
Field_a | Field_b (that will contain random values)
1       | 0.001
2       | 0.005
3       | 0.008

That correlation doesn’t need to be high, I don’t want totally random values. In other words, when Field_a values are increasing, so does Field_b.
Also, both fields may have different value ranges, for example in Field_a values vary between 1 to 10 and in Field_b between 0.1 to 0.5.


texts: How to make the arrows in d.rast.arrow GRASS script visible?

texts:
I'm trying to create a Flow or directional raster map  (with arrows) from an SRTM elevation map that I downloaded from CGIAR.

I already created an aspectMap to be used for the d.rast.arrow. I used the script:
d.rast.arrow map=aspectMap type=grass

and it only displayed a map that is solid green without any arrows.

I also tried the script:
d.rast.arrow map=aspectMap type=grass skip=50

and it displayed dots on the screen for every 50th cell. I already zoomed in to its max, but still I cannot see the arrows.

Is there a problem with my script? or is there an alternative way to create the arrows?


texts: I'm transforming a WGS84 CRS to ETRS 89 using QGis. When I reproject the raster using via raster transformation or via GDAL Warp my pixel change from 30x30m to 26x26m which shouldn't be a problem as far as I know. When I instead save the WGS84 to ETRS89 the pixel get kind of distorted and aren't square anymore but 30x20m. The whole layer looks to be the same size either way. Can somebody explain me why this difference occurs and tell me which transformation is recommended?


texts: Layer definitions for Web app builder layer list widget

texts:
Built a layer list widget using web-app builder and integrated it with our Spring application.

I'm trying to pass state parameter to layer definitions through the following code:

this.operLayerInfos._operLayer.layerObject.layerDefinitions=["STATE IN ('CA')"];
this.operLayerInfos._operLayer.layerObject.supportDynamicLayers = true;
this.operLayerInfos._operLayer.layerObject.disableClientCaching = true;

This is not reflecting in the map.


texts: installing arc.js library in OpenLayers not working

texts:
I have project developed in Angular 5 with Node.js .
Now in my project I want to have flight animation connections and used the example of OpenLayers called 'Flight Animation'. The problem is that there is a library called 'arc.js' used in example which I could not import it. I installed the library using 'npm install --save arc' but still 'arc' is not recognized. My code is :

lineAnimation(){

var style = new Style({
  stroke: new Stroke({
    color: '#EAE911',
    width: 2
  })
});

var flightsSource;
var addLater = function (feature, timeout) {
  window.setTimeout(function () {
    feature.set('start', new Date().getTime());
    flightsSource.addFeature(feature);
  }, timeout);
};
var pointsPerMs = 0.1;
var animateFlights = function (event) {
  var vectorContext = event.vectorContext;
  var frameState = event.frameState;
  vectorContext.setStyle(style);
  var features = flightsSource.getFeatures();
  for (var i = 0; i < features.length; i++) {
    var feature = features[i];
    if (!feature.get('finished')) {
      // only draw the lines for which the animation has not finished yet
      var coords = feature.getGeometry().getCoordinates();
      var elapsedTime = frameState.time - feature.get('start');
      var elapsedPoints = elapsedTime * pointsPerMs;

      if (elapsedPoints >= coords.length) {
        feature.set('finished', true);
      }
      var maxIndex = Math.min(elapsedPoints, coords.length);
      var currentLine = new LineString(coords.slice(0, maxIndex));

      // directly draw the line with the vector context
      vectorContext.drawGeometry(currentLine);
    }
  }
  // tell OpenLayers to continue the animation
  this.Map.render();
};
flightsSource = new SourceVector({
  wrapX: false,
  attributions: 'Flight data by ' +
    'OpenFlights,',
  loader: function () {
    var url = 'https://openlayers.org/en/v4.6.5/examples/data/openflights/flights.json';
    fetch(url).then(function (response) {
      return response.json();
    }).then(function (json) {
      var flightsData = json.flights;
      for (var i = 0; i < flightsData.length; i++) {
        var flight = flightsData[i];
        var from = flight[0];
        var to = flight[1];
        // create an arc circle between the two locations
        var arcGenerator = new arc.GreatCircle(
          { x: from[1], y: from[0] },
          { x: to[1], y: to[0] });

        var arcLine = arcGenerator.Arc(100, { offset: 10 });
        if (arcLine.geometries.length === 1) {
          var line = new LineString(arcLine.geometries[0].coords);
          line.transform(Proj.get('EPSG:4326'), Proj.get('EPSG:3857'));

          var feature = new Feature({
            geometry: line,
            finished: false
          });
          // add the feature with a delay so that the animation
          // for all features does not start at the same time
          addLater(feature, i * 50);
        }
      }
      this.Map.on('postcompose', animateFlights);
    });
  }
});
var flightsLayer = new LayerVector({
  source: flightsSource,
  style: function (feature) {
    // if the animation is still active for a feature, do not
    // render the feature with the layer style
    if (feature.get('finished')) {
      return style;
    } else {
      return null;
    }
  }
});
this.map.addLayer(flightsLayer);

}

I get the error in 'var arcGenerator = new arc.GreatCircle('


texts: Save or add QgsVectorLayer object feature in a layer that is already created

texts:
After getting the QgsVectorLayer result, I want to save its features in a layer named dimensioning, that is already created. How to do it in PyQGIS?

result = processing.run("qgis:explodelines", {'INPUT':temp_lines,'OUTPUT':'memory:'})["OUTPUT"]


texts: Using Focal Statistics to calculate saddle points from DEM in ArcPy?

texts:
I have a conundrum with Focal Statistics and kernel files.  I am been trying to identify saddle points in DEMs.  I am attempting to rewrite this (thanks whuber) in a python script for 10.1.  I have used focal statistics and kernel files for irregular neighborhoods in the past with no issues.  Now, for some reason, I continuously get the following error:

ERROR 010308: Invalid Irregular neighborhood mask

My script is below:

# Import system modules
import arcpy, os, codecs
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")

# Obtain parameters
dem = arcpy.GetParameterAsText(0)
outFolder = arcpy.GetParameterAsText(1)

arcpy.env.workspace = outFolder
DEM = Raster(dem)
kernel = os.path.join(outFolder, "kernel.txt")

# Compute the neighborhood sums involved in estimating the second partial deriviatives.

arcpy.AddMessage('Attempting FocalStat 1...')
n1 = NbrRectangle(3, 3, "CELL")
g1 = FocalStatistics(DEM, n1, "SUM", "NODATA")

arcpy.AddMessage('Attempting FocalStat 2...')
n2 = NbrRectangle(1 , 3, "CELL")
gX = FocalStatistics(DEM, n2, "SUM", "NODATA")

arcpy.AddMessage('Attempting FocalStat 3...')
n3 = NbrRectangle(3, 1, "CELL")
gY = FocalStatistics(DEM, n3, "SUM", "NODATA")

arcpy.AddMessage('Attempting FocalStat 4...')
text = codecs.open(kernel, "w", encoding='ascii')
text.write("3 3\n1 0 0\n0 0 0\n0 0 1")
text.close()
n4 = NbrIrregular(kernel)
gM = FocalStatistics(DEM, n4, "SUM", "NODATA")

arcpy.AddMessage('Attempting FocalStat 5...')
text = codecs.open(kernel, "w", encoding='ascii')
text.write("3 3\n0 0 1\n0 0 0\n1 0 0")
text.close()
n5 = NbrIrregular(kernel)
gP = FocalStatistics(DEM, n5, "SUM", "NODATA")

# Delete kernel txt
os.remove(kernel)

# Compute the determinant of the Jacobian.
# It is multiplied by 10^6 to bring it into a range that SA can display.
arcpy.AddMessage('Calculating determinate of the Jacobian...')
xCell = arcpy.GetRasterProperties_management(DEM, "CELLSIZEX").getOutput(0)
g = (g1 - (gX * 3)) * (g1 - (gY * 3)) / 9.0 - (((gP - gM) / 4.0)**2) / (int(xCell)**4 / 1000000)
g.save("JacobDet.img")

I have tried using the focal stat tool in arc with a standalone txt file with the same results.

Is there something wrong with my code/kernel file?


texts: Getting error in running ExportReport Python code in ArcMap

texts:
I am unable to generate pdf report using ExportReport function of ArcGIS 10.2. getting error of IOError: Could not open report template

I tried different template files but still getting this error.
Template directory is correct. checked all the directory positions.

Code

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Users\Abrar ahmad\Documents\ArcGIS\New_Rwp_Cencus(12-17-2014).mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, "",df):
if lyr.name == "New_Districts_5_without_limit":
    arcpy.mapping.ExportReport(lyr,r"C:\Users\Abrar ahmad\Documents\ArcGIS\test.rlf",r"C:\Users\Abrar ahmad\Documents\ArcGIS\ProjectReport2.pdf","USE_RLF")
del mxd

Complete Error Data:

Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\mapping.py", line 515, in ExportReport
    return report_source._arc_object.ExportReport(*gp_fixargs((report_layout_file, output_file, dataset_option, report_title, starting_page_number, page_range, report_definition_query, extent, field_map), True))
IOError: Could not open report template


texts: How can I calculate the canvas scale in QGIS print composer?

texts:
Using the atlas functionality in QGIS 2.14 The canvas scale varies according to the object in question (national parks).

I want to have a numbered scale based on the canvas print size.  I tried using scale in a text box:

[% $scale %]

This gives me "0" and is not useful. Either this is an error, or I am not using it right.

Requesting pagewidth and dpi works fine. Is there a different way to do this? Based on pagewidth and real world maxx and minx I should be able to calculate the scale. But I can not find any variables for maxx and minx.

How can I present scale on the atlas generated maps?


texts: My code compiles sometime or it gives compile timeout sometime in Google Earth Engine

texts:
I made a project for finding build-up area of a city. It was running fine a month ago but now it compile sometime or give the compile timeout error. Here is my code:

var cloud_masks = require('users/fitoprincipe/geetools:cloud_masks');
var maskClouds = cloud_masks.landsatTOA();

//Loading India image, the extracting data for Haryana (a state in India) and then subsequently Ambala (a district in Haryana)
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11'];
var india = ee.FeatureCollection('ft:1UDdgOCf8DoRJ9bVm-UVbR6CqxtkJToLQjTFd0r0Z','geometry')
.filter(ee.Filter.eq('Name','India'))
.geometry();
var india_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(india)
.filterDate('2014-02-01','2014-09-01')
.sort('CLOUD_COVER')
.map(maskClouds)
.median();
var district = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
.filter(ee.Filter.eq('name','Jodhpur'));
var haryana_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate('2014-02-01','2014-04-01')
.sort('CLOUD_COVER')
.map(maskClouds)
.median();
var input = haryana_image;
input = addBands(input.select(bands));
india_image = addBands(india_image);

//Loading the points from the fusion table and training the classifier
var ft = ee.FeatureCollection('ft:1fWY4IyYiV-BA5HsAKi2V9LdoQgsbFtKK2BoQiHb0');
var ft_builtup = ft.filter(ee.Filter.eq('class',1)).limit(1200);
var ft_nonbuiltup = ft.filter(ee.Filter.eq('class',2)).limit(1800);
ft = ft_builtup.merge(ft_nonbuiltup);
var new_bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11','NDBI','NDVI'];

function addBands(image){
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  var ndbi = image.normalizedDifference(['B5', 'B4']).rename('NDBI');
  var ndwi = image.normalizedDifference(['B6', 'B6']).rename('NDWI');
  return image.addBands(ndvi).addBands(ndbi).addBands(ndwi);
}

// Load a Landsat 8 image to be used for prediction.
var training = india_image.sampleRegions(ft,['class'],30);
var trained = ee.Classifier.cart().train(training, 'class', new_bands);

input = input.clip(district);
input = input.classify(trained);
input = input.expression('LC==1?1:0',{'LC':input.select('classification')});
Map.addLayer(input.clip(district));

/*Export.image.toDrive({
  image: input.clip(district),
  description: 'Test',
  scale: 30
});*/``


texts: Creating geoprocessing package (*.gpk) that contains more than one model?

texts:
I'm making gpk (= geoprocessing package) packages for ArcGIS Runtime 10.2.2.

I'm making several models, and wanted to know if there is a way for including several models in one gpk file or do i need to make several gpks and move between them with code?


texts: Constructing simple zip code polygon map in CartoDB?

texts:
I'm new to using CartoDB.

I want to produce a map that depicts multiple layers, the first of which is simple zip code polygons (i.e. polygons outlining an area defined by a postal zip code).

I would ideally like to stick to using CartoDB to do this.

I then plan to add other layers but need to start with zip codes.


texts: How to get the amplitude in X and Y of an image in Google Earth Engine?

texts:
I have an rectangular image in Google Earth Engine and I would like to know its height and width. How can I do that?

I tried this below, but the result was 0.0013635360788981643
0.5018710804225908 for X and Y, but the amplitude should be something bigger than 0, because I want to divide this image later into a bunch of images with shape 256x256

listaCoord = ee.List(geometria.coordinates().get(0))
t0 = listaCoord.get(1)

t1 = listaCoord.get(3)

x0 = t0.getInfo()[0]
y0 = t0.getInfo()[1]
x1 = t1.getInfo()[0]
y1 = t1.getInfo()[1]

ampX = (x1 - x0)

ampY = (y0 - y1)


texts: Can a WMS layer be added as a feature layer to CartoDB?

texts:
There are plenty of pointers that show how to add WMS as a base layer in CartoDB. However, I want to add a WMS layer to the [Layer Switcher] so that it is available as a queryabl feature layer.

I was hoping this option would be available through "Add new Table" and for a source of the table we would see a WMS connector. I see ArcServer is available as a layer provider, but no generic WMS.

Am I missing something or is it simply not possible to add a WMS layer directly to a CartoDB visualisation? I could not find the point covered in their docos...

The implication of course is that we will need to somehow script an update, export/import pipeline to get some of our (often updated) data out of Geoserver and into CartoDB. Seems like a shame when there's a lovely WMS endpoint just sitting there ready to be exercised by CartoDB...if only it were possible??


texts: How to pass polygon to query.geometry with ArcGIS API for JavaScript

texts:
I have to highlight district boundary and corresponding point. So I am trying to do like below code. But I don't know how to pass district polygon geometry to query.geometry.

var queryTask = new QueryTask(".../arcgisserver/rest/services/CRD/CRD2/MapServer/0");

var query = new Query();

query.geometry = event.mapPoint; //your district polygon geometry
query.outFields = ["*"];
query.returnGeometry = true;

//Execute task
queryTask.execute(query, showResults);

function showResults (results) {

    var resultItems = [];
    var resultCount = results.features.length;

    for (var i = 0; i < resultCount; i++) {
        var featureAttributes = results.features[i].attributes;
    }

}

Is it possible to pass polygon geometry into query task. If so how can do it.


texts: Loading Historical Google Satellite Imagery with QGIS QuickMapServices plugin?

texts:
Loading older google satellite image with OpenLayers plugin in QGIS? is a similar question that was asked 4 years ago.

Is QGIS QuickMapServices able to load Google Earth's (aka Google Satellite's) historical imagery?


texts: Publishing commercial ArcGIS Server Geoprocessing Services?

texts:
I've developed a set of ArcToolbox Python scripts, which I've been using on certain projects for clients. The approach we've been taking is to have our clients pay us to run these tools, and provide data and a report at the end of the project.

This approach is working well so far, but has obvious limitations. I've been thinking that there could be advantages to distributing the tools - we'd be able to reach a wider market, and it could drive more consulting work. I think best way to do this might be to publish the toolbox to ArcServer and charge clients for access (via ArcMap), but I'm hardly familiar with ArcServer, so I wanted to get clarification on some things.

Is it actually legal to charge users for access to ArcServer geoprocessing
services via ArcMap? I don't see how this would be different than charging for access to a web map with geoprocessing functions, but I wouldn't want to violate ESRI's terms of service.
When a Python script is published on ArcGIS Server, can ArcMap clients view the code, or is it hidden? If we publish this set of tools at all, I don't want users to be able to view or copy the code.
Is this even a good idea?

Original post on ArcGIS forums: http://forums.arcgis.com/threads/62752-Geoprocessing-Service-General-Questions


texts: How to keep track of changes to geo-referenced data and maps?

texts:
At my company we are starting to get more serious about keeping track of changes to our geo-referenced data and maps.  I'm looking into this and am wondering how people do it at other companies.

We have decades of historical data that is often suspect, and gets geo-referenced multiple times as new information comes to light.  We want a document that is used for each transformation to record the source (coordiantes, topography, other map), projection, etc, etc.  Our GIS team is 5 people and growing with 50 in house people accessing the data as well as outside organizations and government.  At the moment we use Windows 7, ArcGIS Desktop 10.0, Geosoft Montaj, AutoCAD 2012, Microsoft Access '07, and Sharepoint, as our primary programs to make changes.

Is there an automatic logging program that records everything that was done during a session in Arc?  Which would also keep track of projections and transformations, point creation, etc.  Or is this thing best done in Notepad by hand?

Thanks for your feedback.


texts:
I have been successfully using the Gdaltools module in Python to export from PostGIS to shapefile using:

import gdaltools

ogr = gdaltools.ogr2ogr()

conn = gdaltools.PgConnectionString(host=HOST, port=5432, dbname=db, schema=scema, user=user, password=PASSWORD)

ogr.set_input(conn, table_name="tablename", srs="EPSG:4326")

ogr.set_output('test.shp')

ogr.execute()

I would like to add an SQL query to the mix (since the whole table is currently exporting) and not sure where it should go.


texts: I am running a Euclidean Distance analysis on a large file. I have a polygon layer that includes all Urban areas in Canada, US, and Mexico merged into one multipart polygon file. I am attempting to calculate Euclidean Distance at the 10 x 10 m pixel resolution. It has been three days, the hour glass in the results window is still there and the status bar at the bottom is still running with "Euclidean Distance...". I also set a processing extent to just outside the coastal boundary of the entire North America. Is there a way to see how much progress has been made, or how much time (or percent) is left to do? i have no idea if this is going to run for another hour or another 3 days.


texts:
I try to use qgis to create a custom topographic map for outdoor use: UTM, 1:50000.

What these maps usually have is a 1km x 1km grid, which at the bottom is labeled with a shorthand for the corresponding UTM grid coordinate:

Say the UTM easting of a vertical line is 05532000, then I'd label the line with 32. The 055 is redundant because every line in the map will have that part (provided the map doesn't cover too large an area), and the 000 is redundant as the grid lines are placed at multiples of 1000.

Now, I used qgis to create the grid, and I added a column to the corresponding data file containing the labels I'm after, but when I use the "Display labels" functionality in the properties of the shapefile, the labels occur in the center lines of the map, whereas I'd like to have them at the margin. I didn't really find a way in qgis to do that and would appreciate any help.

(If other open source software is better suited to create nice topographic maps from a GeoTIFF file, I'd also appreciate that, as long as it runs on Linux or Mac)


texts:
Issue : AddIn button/tool force Arcgis 10.2.2 to close.

Development environment: Visual Studio 2013, .Net framework 4.5 (I tried with 4.0 also), Arcgis for desktop 10.3
Works fine on development machine with version in AddIn config file as 10.1 or 10.2 or 10.3

User environment: Arcgis 10.2.2. AddIn get installed with version in AddIn config as 10.1 or 10.2. I am able to see the tool bar and tools. The moment I click in map view after taking a command from my AddIn Arcmap window exits without any message


texts:
I have created a grid, which I want to use in an ArcGIS Online application. To create the grid (fishnet) I have used the local projection, which is in meters. I avoided using WGS84 as it is in degrees and is not so accurate. As you can see in the attachment the grid looks ok but as expected the squares are not in the same line. Since I will use this map for reporting purposes, the tilted squares don't look so good.
Do you have any ideas how I could "fix" this in order to look more neat?


texts: QGIS filling more than one column in data attribute table at once

texts:
I would like to fill up the columns, which values are the same. There are 2 columns (in future it can be more).
Unfortunately, I've encountered a problem, described in the image below. I don't know really how to solve it.

I've marked the stuff I can put in. The column data should be integer here. I don't understand why it doesn't work.
I found the similar problem under the thread below:
QGIS Field Calculator Error for String values
which applies to strings. Since I quote my number it, the quote disappears and problem remains the same.
Is there any way to use the field calculator for more than one column at once?


texts:
I have issues in re-projecting a vector layer to match the rest of my data. I have a polygon vector layer with districts, as well as a point layer in WGS84. I know these are correct since I've been viewing them correctly in Google Maps and Google Earth in KML format.  Then I have a vector layer of lines (given to me by another party) that appear far away from my other data in QGIS, and also in a completely different scale (the lines are way larger in comparison to by districts). The layer spatial reference units are also different in this layer, it is not shown in degrees as it is in the rest of my data.

I have been trying to re-project this to WGS84 (EPSG:4326), but the new layer appears in the same spot as before. I have also tried re-projecting it to WGS84 Pseudo Mercator (EPSG:3857), when I do that, the data disappears (can't find it though I try to zoom to it).

How can I alter/re-project the line-layer to match the rest of my data?


texts: Raster layer with unique vector for each raster cell?

texts:
I have been looking around several hours, even days, but couldn't find a clear answer to my question.

I am working on a WebGIS tool that interprets local geology. From more than 70 asc-files, I would like to create one raster layer that contains in each raster cell the local geology, e.g. as a 2 x n vector with a code for the present formations and their thickness.

First of all, is this a good idea? Or is it better to query all layers in the WebGIS tool for the specific location?

I tried already to create a raster file, where each cell value was linked to a vector in a sort of attribute table. However, I was unable to save this table in Python, unless as a txt-file, and directly link it to the raster file. Is there a specific raster type that does this?

Thank you!


texts: Bookmarking the rotation, scale, and extent in QGIS

texts:
I have one map that I want to bookmark at multiple
rotations/orientations.

However, when I use the bookmark tool to capture the view, I am only able to capture the scale and extent of each view, and not the rotation.

For example, I saved this view as a bookmark, labelled 'Horizontal" (note the rotation at -29.0):

When I change the rotation to 60 degrees to get the following view:

And I use aforementioned 'Horizontal' bookmark to get back to the original view, the scale and extent of map changes, but the rotation remains at 60, instead of changing to -29:

I have tested this on QGIS versions 3.0, 2.99, and 2.18.16 on Windows computer.

Is there a better way to create multiple default layouts that include map extent, scale, and rotation?


texts: Same confusion matrix depiste hyperparameter tuning using TrainImageClassifier in OTB

texts:
I am trying to run a land use land cover classification using TrainImageClassifier in OTB. Here are the settings I used :

The problem is that no matter how much I tweak the hyperparameters of any of the available algorithms I end up with the same final confusion matrix and scores:


texts:
I am new to python and QGIS (have some experience with arcpy). I am trying to multiply two .tif files together using gdal_calc.py but I can't get it to work. It just returns the error invalid syntax. Here is the code. I am doing this in the QGIS Python command line.

gdal_calc.py -A people.tif -B places.tif --outfile=c:/test.tif --calc="A*B" --NoDataValue=0

I then tried this after reading an answer to another question, which ran and printed a number. However, it has not output a .tif file:

import os

gdal_calc =  'python C:/wheregdal_calc/islocated -Apeople.tif B-places.tif --outfile=c:/test.tif --calc="A*B"'
os.system(gdal_calc)


texts: Is there a limit to the amount of data in a spatial join in QGIS?

texts:
I am using Quantum GIS and trying to spatially join a point grid with polygons that are overlapping. If I get the sum of intersecting features, then I can get a good raster value. Well, every time I attempt the spatial joing with the summary, I get the following python error:

  An error has occured while executing Python code:

  Traceback (most recent call last):   File
  "/usr/share/qgis/python/plugins/fTools/tools/doSpatialJoin.py", line
  99, in accept
      res = self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar)   File
  "/usr/share/qgis/python/plugins/fTools/tools/doSpatialJoin.py", line
  214, in compute
      if k == "SUM": atMap.append(sum(numFields[j])) TypeError: unsupported operand type(s) for +: 'int' and 'QPyNullVariant'

  Python version:
  2.7.4 (default, Sep 26 2013, 03:23:21)  [GCC 4.7.3

Any idea why I may be getting this error? I am trying to join 65,527 points to 121 polygons in this case.

Thanks for the help!


texts: Importing metadata from HTML into ArcCatalog fails

texts:
I downloaded restricted areas shapefile from here:

https://webgis.dme.qld.gov.au/webgis/webqmin/viewer.htm

It came with a rst.html file, which, on opening, displayed a page exactly the same as the 'Description' tab in ArcCatalog. Enthused, I imported this file into ArcCatalog, hoping it would transfer all the information automatically. Unfortunately, however, it only imported the bounding box coordinates and put the name of the shapefile as the heading.

Anyone know what went wrong?

Wil


texts: Definitely delete more than 1000 images from an image collection in Google Earth Engine using the code editor!

texts:
I need to definitely delete more than 1000 images from an image collection in Google Earth Engine using the code editor.
After some attempts I arrived to the following code:
var deleteA = function(image) {

    return ee.data.deleteAsset(image,true);

};

print(collection.map(deleteA));

But it didn´t work. I dont´t understand the reason. Do you what's wrong in my code?


texts: How to access the no data value of a raster with GDAL?

texts:
Is there a similar way as in ArcPy

import arcpy

desc=arcpy.Describe('rasterName')
print desc.noDataValue

to access the no data value using GDAL?


texts: QGIS select points within certain radius

texts:
I've seen other threads ask very similar questions, but I haven't seen one that answers this specifically:

So, I have 4 Total sets of points.  1 set is a set of addresses in Arkansas, the other 3 are different types of toxic waste contamination sites.  What I am trying to find out is how to run a query to find and select which addresses are within 20 kilometers of any of those 3 contamination sites.

Does anyone know how I could go about doing that?


texts: How to include an overview map control code in qgis2web code

texts:
Please, forum I would like to know how to add an overview map control code in my qgis2web  generated OpenLayers code.If am not clear enough the link below demonstrates what I mean
https://openlayers.org/en/latest/examples/overviewmap.html


texts: Georeferencing DXF file to use in Google Earth

texts:
I have read up on several questions that seen to have similar problems in georeferencing dxf files from CAD. None of them are clear, so I am trying here.

I have CAD files of solar plant layouts, roads, fences, etc in DWG format. I have used Any DWG DXF Converter to convert into DXF, which I can open in QGIS. The problem is there are no recognisable scales or coordinate system. The scale is currently 1:377579882 when viewing the plant layout from just above and the coordinate on the mouse point is 95675, -3008800.

I have tried reading about the Affine plug in, but can't work out how to use it as the the Transformation Matrix makes little sense without explanation. I have tried using and reading up on the V.Transform module in GRASS, but I don't know what the "Shifting Value for X and Y coordinates" is or the "Scaling Factor for X and Y coordinates" - None of the previous answers seem to actually tell users how to run these modules or find the correct inputs.

The purpose of this exercise is to take the CAD files (which I can make into DXF), load them in QGIS (possible convert to Shapefile if necessary) and save them as KML files to open in Google Maps.


texts: How to calculate total area for classified polygons?

texts:
I am trying to calculate the area of landuse types I have classified into 4 classes. I currently have the area for each polygon in the Shape_Area field and the respective class in the gridcode field.
I'd like to know what steps can be done to sum the areas for each class, as I've forgotten how its done.


texts: Why are the output point shapefile PRJ files empty?

texts:
I have tried everything I could find (using C++)  to create a new .prj file for an output point shapefile using OGR.  While I get no compilation or execution errors and ArcMap accepts the file as a point file, ArcMap still claims that my output shapefile has no spatial reference:

Data Type:  Shapefile Feature Class
Shapefile:  G:\GIS\Working\FIrani\Ubuntu\0_From_Ubuntu\Seeds001.shp
Geometry Type:  Point
Coordinates have Z values:  No
Coordinates have measures:  No

Coordinate System:  <Undefined>

How do I get OGR to actually apply my ACEA specified coordinate information to its output shapefile?  Do I actually have to delete the current .prj and create and write out a new .prj file myself?

My current (kitchen sink) approach is as follows:

    /*
    Define projection for output seed file.
    --------------------------------------- */
    char *pszWkt = NULL;
    oSRS.SetProjection("ACEA");
    oSRS.SetProjCS("ACEA");
    oSRS.SetWellKnownGeogCS("ACEA");
    oSRS.SetACEA(29.5, 45.5, 23.0, -96.0, 0.0, 0.0);
    oSRS.morphToESRI();
    oSRS.importFromWkt( &pszWkt );
  pioSeedLayer = pioSeedDS->CreateLayer( pszSeedsfn, &oSRS, wkbPoint, NULL );
  // Ouput projection via oSRS.
  if( pioSeedLayer == NULL )
  {
    printf( "Layer creation for seeds failed.\n" );
    exit( 1 );
  }


texts: Splitting single raster within multiple polygons using ArcMap

texts:
I have a raster ('slope') which content only one attribute in the table. On the other hand I have a feature which content several polygons (catchments) in the attributes table. I would like to obtain the area or count of the 'slope' for each catchments.

Is there any faster way to do it without separating the catchments and masking the rasters one by one?


texts:
When adding an image overlay (specifically a GeoTiff exported from QGIS), the image looks great until you turn on 3D buildings or terrain. These features pop through the image overlay, and disrupt the image. In an ideal world, I would want the image to wrap around the 3D features.

Use Case: With an elevation map, it would be helpful to see 3D elevation features coinciding with the an image overlay colored to indicate where the high elevation is.


texts: Applying Conditional Symbology in ArcMap?

texts:
I'm in charge of a database that contains property information.  I would like to symbolize parcels of land that we have a 'written agreement' with the owner differently than parcels that we do not have written agreements with.

My solution was to create a calculated field (1 or 0) based on if there was any data in the "Agreements" field.

  I'm looking for a more dynamic solution (not having to re-run the
  field calculation) so that when an agreement is entered into the db,
  the symbology automatically changes.

This is important because our field staff take the db out into the field on hand-held devices; they need to know which properties they are allowed to access (green) vs. which properties they cannot access (red).


texts: Adding tracks after importing GPS time/position into QGIS?

texts:
In properties there is "Joins".

I was wondering if this would give me lines linking the points (time order) but can't see how?


texts: FieldMap in Arcpy Merge for Selected Fields Only

texts:
I am trying to use FieldMap in arcpy to merge only a few of the same fields from 8 feature classes. There are many different fields in all of the different feature classes, but I have identified the 9 fields I want in my final merged feature class. Instead of merging all the fields I identify it only chooses the first field in my list and I end up with a merged file with one field "COUNTY_F". What am I missing?

fClasses = []
i=0

#Loop through MyGDB.gdb and add each feature class name to array
fcs = arcpy.ListFeatureClasses()

for fc in fcs:
    fClasses.append(fc)
    i = i+1

print fClasses
#Shows as "Clallam_WorkingForest_2013, GraysHarbor_WorkingForest_2013, Jefferson_WorkingForest_2013, King_WorkingForest_2013, Kittitas_WorkingForest_2013, Mason_WorkingForest_2013, Pierce_WorkingForest_2013, Snohomish_WorkingForest_2013"

#Map Fields for Merge
inFields= ["COUNTY_F", "UNIQUE_ID", "TAXPAYER", "PARCEL_ID", "CURRENTUSE_CODE", "CURRENTUSE_DESC", "PRESENTUSE_CODE", "PRESENTUSE_DESC", "ACRES_ASSES"]

fm = arcpy.FieldMappings()
fm_type = arcpy.FieldMap()

for field in inFields:
    fm_type.addInputField( "Clallam_WorkingForest_2013", field.name)
    fm_type.addInputField( "GraysHarbor_WorkingForest_2013", field.name)
    fm_type.addInputField( "Jefferson_WorkingForest_2013", field)
    fm_type.addInputField( "King_WorkingForest_2013", field)
    fm_type.addInputField( "Kittitas_WorkingForest_2013", field)
    fm_type.addInputField( "Mason_WorkingForest_2013", field)
    fm_type.addInputField( "Pierce_WorkingForest_2013", field)
    fm_type.addInputField( "Snohomish_WorkingForest_2013", field)

fm.addFieldMap(fm_type)
print "Fields mapped!"

arcpy.Merge_management(fClasses, "Merged", fm)
print "Merge: SUCCESS."


texts: Samples of US Residential or business addresses without zip+4?

texts:
I'm looking for a few deliverable residential or business addresses without a zip+4.  By "residential or business", I mean, not a PO box and not military addresses.

According to this answer on Ask.com (which no longer exists), such addresses do exist, but I'd like to see some samples.


texts: Convert PostgreSQL Polygon to PostGIS POLYGON

texts:
I have an existing table in postgres which contains a plain old posgresql polygon.  I have added a new column using

SELECT AddGeometryColumn ('public','fences','geom',4326,'POLYGON',2);

I would like to populate this GIS column with the polygon data from the existing column.

The existing column is called poly and the data looks like the following:

"((-74.026142,40.636264),(-74.026142,40.636438),(-74.025818,40.636438),(-74.025818,40.636264),(-74.026142,40.636264))"

Is there a query I can run that will pull that poly field convert it to a postgis polygon and insert it into the new geom column?


texts: I heard today that the Arctic Circle, along with it's siblings the Antarctic Circle and the Tropics of Cancer and Capricorn, are different from year to year. Is that true? Can you please explain how/why this is so? I'd always thought they were defined by particular latitudinal parallels. e.g. 66°N or S etc.


texts:
I'm trying to use Clip in QGIS with a polygon layer (buffer of 5000 meters) as input and another layer as clip. However, the results are just empty with no attribute table and no image results.

I've succeeded with lesser buffers (20, 50, 200 meters etc.) but it seems like when I have a buffer that's quite large it doesn't work. The buffers are from about 18 000 points so it's heavy data too.

They have the same reference system and all, it worked as I said with lesser buffers but not for this one.

ArcMap error says "000117 : Warning empty output generated."

I also tried this in ArcMap but that just generates and error. I can't wrap my head around this, so hopefully someone can have some input.


texts: Using formatting option 'aligned' in 'to_dms' function in QGIS 3.4.2

texts:
Using formatting option 'aligned' in the expression:

to_dms(x(transform($geometry,'EPSG:32632','EPSG:4230')), 'x', 2, 'aligned')

not only is the alignment of degrees, minutes, seconds as expected, but also the suffix E (which should be added with the 'suffix' format option) is added.

Is there a way to eliminate the suffix?


texts: Google Earth Engine: Cannot read property 'percentile' of undefined

texts:
I have a problem running the code below. I'm using the API through Node js. I receive

  TypeError: Cannot read property 'percentile' of undefined

I don't understand why. Can you help me, please?

function test() {
      ee.data.authenticateViaPrivateKey(config.gee, async () => {

          console.log('init')
          var geometry = await ee.Geometry.Polygon([[[-58.22065005792236, -38.39475941140011],[-58.21198115838623, -38.391698569574615],[-58.222109179626464, -38.38813303001853]]])

          var collection = await ee.ImageCollection('COPERNICUS/S2').filterBounds(geometry).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 5));

          var cloudMasked = await collection.map((image) => {
            var cloud = ee.Number(2).pow(10).int();
            var cirrus = ee.Number(2).pow(11).int();
            var qa = image.select('QA60');
            var mask = qa.bitwiseAnd(cloud).eq(0).and(qa.bitwiseAnd(cirrus).eq(0));
            return image.updateMask(mask);
          });

          var ndvi = await cloudMasked.map((image) => {
            return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
          });

          var reducer = ee.Reducer.percentile([90]);
          var ndvi90 = await ndvi.select(['NDVI']).reduce(reducer).clip(geometry);

          var url = await ndvi90.getDownloadURL({'scale': 10})

          console.log(url);

      });
}

Is the use of async/await correct?


texts: Edit raster classification manually using polygons in ArcMap?

texts:
I want to edit the elevation/z values of certain zones on a raster in ArcMap. I followed this example by creating a polygon, assigning a value of 7.5, then converting the polygon to a raster, finally using the raster calculator to create a new layer:
Manually edit a raster to improve a classification?

However, after using the raster calculator I was left with a raster containing the values of the polygons only, with the remaining values from my original raster as 'nodata'.


texts: Trying to repair layer paths on opening of a QGIS project from Geopackage

texts:
It's some time ago, that I was trying to use the bad layer handler in QGIS to change some orphaned layers during project opening. Unfortunately, I was not able to call the default bad layer handler in advance, to let the user manually handle all the remaining unhandled layers. Honestly, I got stucked.

Yesterday I accidentally found QgsPathResolver, which looked very promising to me. I was very delighted as I saw the examples in the API description.

Some background: my main aim is to automatically repair a QGIS template project, that was stored together with all it's layers inside of a GPKG file. Everytime the GPKG file is copied and renamed, the relative layer paths are not automatically updated.

Unfortunately, I couldn't achieve a replacement of the old GPKG filename inside the relative paths (see code below).

Can someone shed some light on this? Do I face a bug here or a missing implementation?

BTW: I'm using QGIS 3.12.0 (Windows 10)

from qgis.core import QgsPathResolver,QgsProject
import re

def my_pathpreprocessor(path):
    projectFileName = QgsProject.instance().fileName()
    newPath = path
    if '.gpkg|' in path and path.startswith('./'):
        gpkgProject = re.findall(r'[^/\\]+\.gpkg',projectFileName)[0]
        gpkgProjectPath = re.sub(r'^[^\:]+\:(.+/)[^/]+\.gpkg.+',r'\1',projectFileName)
        print('Old path: '+path)
        if not gpkgProject in path:
            newPath = re.sub(r'[^/\\]+\.gpkg',gpkgProject,path)
            newPath = gpkgProjectPath + newPath[2:]
            print('New path: '+newPath)

    return newPath

QgsPathResolver.setPathPreprocessor(my_pathpreprocessor)


texts: Why are .grd files used? (Mining and construction industries)

texts:
I'm researching the .grd file format and its not clear to me why they are used, especially regarding the mining and construction industries. I understand they demonstrate changes in elevation or other geographic attributes in a binary format.

What analytical operations are .grd files used for in the mining and construction industry?


texts: FME Server not able to use sde files after ArcGIS Server Update

texts:
This is the description of the environment and situation:

Environment:
Windows Server 2016, 64 Bits
ArcGIS Products:
ArcGIS Web Adaptor 10.7.1,
ArcGis Server 10.7.1,
Python Additional packages:
Pycriptodome 3.9  https://pypi.org/project/pycryptodome/3.9.8/
FME Server:
FME Server 2019.2.3.2, Express install
Build 19825 - win64
2 local Engines, running under the Local System account.

There are 3 workspaces that include the reading from a geodatabase. All workspaces have been working without any issue.
The ArcGIS products were updated to 10.8.1 (Server/Web adaptor) following the instructions of ESRI.
After that procedure the workspaces are failing showing an error accessing to an sde file created  during the execution of the workspace.

The log in the 3 workspaces is similar.
The workspaces were created with FME Workbench 2018.1.
In order to identify the problem, FME Workbench 2019 was used to create a new workspace that uses one of the generated sde files to connect to your SDE Geodatabase: It works with out any issue.
Then, publish that simplified workspace along with the sde file to FME Server and run it there: It works with out any issue.
One piece of one of the workspaces was rebuilded using FME workbench 2019, it includes the generation of the sde file and reading and writing to the database, It works, no issue at all.
Some colleages say, that it is usual that workspaces get corrupted, and you must rebuild them.
The problem here is that those Workspaces were not touched in the update process of the ArcGIS Server, so, how they got corrupted? Moreover, those are huge workspaces, rebuild them will take about 3-4 weeks. (I already try a big copy and paste into a new workspace, It did not work :( ).
Is possible to fix the Workspaces in an acceptable time frame?


texts:
I am trying to use the pysheds package to calculate flow accumulation from a DEM. However, when reading the DEM .tif-file, I get the following error. I can't find any more information as to what this means. My code is as follows:

home_dir = os.path.dirname(os.path.realpath(__file__))
test_file = os.path.join(home_dir, r"N040W010_N045W005\repr-N042W006_AVE_DSM.tif")
grid = Grid.from_raster(test_file, data_name='dem')
plotFigure(grid.dem, 'Elevation (m)')

The error I get is this:

  CRSError: Invalid CRS input: CRS({'init': 'epsg:32630'})

I would love to mean what this means, I have changed the projection from a global to a local one, but with the global projection, I had the same error.


texts: QGIS create line-layer from points, sorted by date/time

texts:
I have point data from telemetrie of a raptor-bird. The point-layer has id, longitude,
lattitude and datetime-fields.
Now I want to create a line-layer which connects all the points, sorted by datetime. So for example the first line would connect the points 1 and 2, the second line the points 2 and 3 the third line the points 3 and 4 and so on.
As addtitional attribute-data each line should contain the datetime-value from its starting point.

I alread tried the "points to path" algorithm, but this gives me either a single polyline (I need seperate lines for each point-connection) or, if I try to group the lines for example by datetime, multiple single lines with the length of 0.
I also used the algorithm explode lines on the polyline/path, which gives me the needed single lines, but I cannot figure out how to give them the datetime-attribute from their starting point.
All referred QGIS-Plugins (points2one, pointConnetr etc) I found on StackExchange doesn't work on QGIS 3.x or are discontinued at all.


texts: Is it possible to maintain primary keys on versioned featureclasses without arcobjects code?

texts:
Since it is not a good idea to use ObjectIDs as primary keys in relationships, some other column will need to be used.  Is it possible to use DBMS sequences to populate primary key columns in versioned featureclasses, or will I need to write arcobjects code to do this?

I've heard Global IDs can be used, but non-GIS DBA's maintaining related tables often frown at this idea.


texts: Intersection throws com.vividsolutions.jts.geom.TopologyException: Found null DirectedEdge

texts:
When I make union of 379202 List<Polygon> rectangles rectangles and execute intersection with bounds UnaryUnionOp.union(rectangles).intersection(bounds) then I get:

com.vividsolutions.jts.geom.TopologyException: Found null DirectedEdge
    at com.vividsolutions.jts.geomgraph.EdgeRing.computePoints(EdgeRing.java:139)
    at com.vividsolutions.jts.geomgraph.EdgeRing.<init>(EdgeRing.java:68)
    at com.vividsolutions.jts.operation.overlay.MaximalEdgeRing.<init>(MaximalEdgeRing.java:65)
    at com.vividsolutions.jts.operation.overlay.PolygonBuilder.buildMaximalEdgeRings(PolygonBuilder.java:104)
    at com.vividsolutions.jts.operation.overlay.PolygonBuilder.add(PolygonBuilder.java:78)
    at com.vividsolutions.jts.operation.overlay.PolygonBuilder.add(PolygonBuilder.java:67)
    at com.vividsolutions.jts.operation.overlay.OverlayOp.computeOverlay(OverlayOp.java:256)
    at com.vividsolutions.jts.operation.overlay.OverlayOp.getResultGeometry(OverlayOp.java:189)
    at com.vividsolutions.jts.operation.overlay.OverlayOp.overlayOp(OverlayOp.java:92)
    at com.vividsolutions.jts.operation.overlay.snap.SnapIfNeededOverlayOp.getResultGeometry(SnapIfNeededOverlayOp.java:96)
    at com.vividsolutions.jts.operation.overlay.snap.SnapIfNeededOverlayOp.overlayOp(SnapIfNeededOverlayOp.java:58)
    at com.vividsolutions.jts.geom.Geometry.intersection(Geometry.java:1338)

Union of rectangles:

Bounds for intersection:

Why do I get the exception for the intersection? How to investigate the problem?


texts: Stratified random sampling with Google Earth Engine Error

texts:
I am trying to sample points from a classified raster using Google Earth Engine, but am running into some problems. I developed a script based on the Classification video located here.

Here is a toy script that attempts to randomly sample a stratified raster.

var ROI = ee.Geometry.Rectangle(-82.562342, 35.589551,-82.562342, 35.589551); // define area

var features = [
    ee.Feature(ee.Geometry.Rectangle(-82.56123, 35.59053,-82.55213, 35.59674),{class: 0}),
    ee.Feature(ee.Geometry.Rectangle(-82.54245, 35.591453, -82.537172, 35.596496),{class:1})
];

var vect = ee.FeatureCollection(features);
print(vect);
Map.addLayer(vect);

var image = vect
  .filter(ee.Filter.neq('class', null))
  .reduceToImage({
    properties: ['class'],
    reducer: ee.Reducer.first()
});

Map.addLayer(image,{min:0,max:1,palette: ['black','white']});

var counts = {
  0:1000,
  1:1000,
};

print(counts)

var stratified = ee.Dictionary(counts)
    .map(function(klass,count) {
    klass = ee.Number.parse(klass)
    var masked = image.updateMask(image.eq(klass))
    return masked.addBands(ee.Image.pixelLonLat())
        .sample({
          region: ROI,
          numPixels: 5000,
          seed: klass})
        .randomColumn('x')
        .sort('x')
        .limit(ee.Number(count).min(100))
        .map(function(f) {
            var location = ee.Geometry.Point([f.get('longitude'),f.get('latitude')])
            return ee.Feature(location,f.toDictionary())
            })
    }).values()

stratified = ee.FeatureCollection(stratified).flatten();

print (stratified.reduceColumns(ee.Reducer.frequencyHistogram(),['Class']));
Map.addLayer(stratified);
print(stratified)

Unfortunately, I receive the error message: The default WGS84 projection is invalid for aggregations. Specify a scale or crs & crs_transform.

This error is mentioned on the projections page, but, doesn't offer any real solutions. Does anyone know how to aggregate using the aerial imagery in WGS84 on Google Earth Engine?


texts: Drawing tools to add circles, arrows, lines for QGIS?

texts:
I'm pretty new in QGIS and i'm searching a drawing tool like it is possible to find in microsoft word or excel.  I search something very user friendly (for dummies!) to add circles, arrows, lines... with no need to create a new layer!!

Can you help me?


texts:
I'm trying to find the shortest path between two points (green and red points) on a road layer using QGIS 2.12.3. The roads layer has appropriate CRS.

Before sign the "start" and "end" points I approximate to the crossing lines, in order to sign exact points

But I get no result


texts: How to export lines from ArcGIS to be used in QGIS

texts:
After a big bunch of bureaucracy we have been giving green light to request government geodata.  The requested data includes mostly seperate landplots, roads and rivers.

We originally were told we'd receive the data in DWG, but we've found out that the government agency is using ArcGIS.

We're using QGIS, and we were going to import/convert the DWG's again - that's double work.  Therefore we'd like to reform our request by asking for ArcGIS data rather than DWG's.

Local bureacracy won't allow us to sit down with the actual people doing the export.  We have to write down our request.

What do we ask them?  What method would be easiest for them to export (ArcGIS) and for us to import (QGIS)?


texts:
Is it possible to assign a keyboard shortcut to a (custom or not...) script of the Processing Toolbox in QGIS?
I have not found any information about it so far.


texts: Viewshed 2. Understanding of VAT and observer-region relationship table

texts:
I am having a hard time using new Viewshed 2.
I am trying to join values from result raster with observer-region relationship table (ORRT) using Value from raster and Region from ORRT.

I couldn't find any workflow of connecting VAT and ORRT. Here is my approach:

a. Convert raster to polygon

b. Inner join ORRT using polygon.Region = VAT.Value

c. Calculate area for polygon

d. Create summarize table to show observers with areas

However, I am a bit concerned of using it since there are a couple of issues you can see below.

Concerns:

1. There are duplicating Values in VAT.

2. There are negative Values in VAT and no corresponding pixels with that Values in the raster.


texts: QGIS project in Lizmap showing only points vector layers

texts:
I am trying to set up a QGIS project in Lizmap, but the problem is that I can't see polygons or lines vector layers.

They appear in the legend; I can zoom on them or use "select by location"; and I can edit them: I can change the attributes and even add new features. The features I add in Lizmap are correctly displayed if I open my layers in QGIS using WFS/WMS or PostGIS, and are saved in my attribute tables. But they are not visible in Lizmap (they "disappear" as soon as I finish drawing them in Lizmap).

The problem doesn't arise with "point" geometry: my points vector layers are correctly displayed in Lizmap, even if they have a complex symbology (heatmap). Raster layers (eg orthophotograph used as base layer) also display correctly.
I did a few tests, and I have the same problem for shape layers used as base layers: they only show if their geometry type is "point".

I have no idea where the problem comes from.


texts: Statistical analyses for survey against raster?

texts:
I am somewhat of a newbie to geostatistical analyses, but pretty
familiar with GIS and remote sensing packages. I was wondering if
anyone had suggestions on appropriate geostatistical analyses for a
particular data set. The data set consists of a large-scale survey at
the land plot level consisting of questions related to political
participation (an ordinal scale), and raster data from a binary
landcover classification of forest cover. The survey data I am envisioning joining to a .shp file of landplot boundaries. I had wondered about either converting the vector to raster and analyzing raster against raster, or converting the landsat landcover map into a vector and analyzing vector x vector.

The research question is whether political participation (potentially
plus other variables) affects the presence of forest cover. What I´m
imagining is some sort of multivariate linear regression in terms of
the survey data, but when it comes to adding in the values from the
raster I´m somewhat stymied. In doing some background research on
potential approaches, it seems like looking at discrete variation
might be the way to go.

I´m using ArcGIS 10 and Erdas 2011.


texts: Populate field in one layer with corresponding value from another layer

texts:
I have 2 existing vector layers: 1 point, 1 polygon. In almost all cases the points overlap a polygon.

The point layer has multiple features (1,000+) with an existing field which I need to retrospectively populate with the equivalent value for the polygon feature which they overlap. How do I do this (without doing it manually)?

In future ideally I would like this to happen automatically when features are created/edited.

I have looked at various similar threads on here, but none of them quite answer the question. It appears that the Autofields plugin might be the solution, but that does not seem to be available anymore?

I am currently using QGIS 3.6.3


texts: What has replaced QGIS SPIT (Shapefile to PostGIS Import Tool)?

texts:
On page 132 of PostGIS in Action, Second Edition it says,

  QGIS SPIT QGIS comes packaged with a tool called SPIT that allows you to batch-load Esri shapefiles into a PostgreSQL database, bypassing rendering. One crippling annoyance with this tool is its automatic conversion of all field names to uppercase. In PostgreSQL, uppercase field names must be double-quoted in queries, and that’s a lot of extra typing of quotes. If you’re going to use SPIT, we suggest that you create a function to convert everything back to lowercase after import.

I found it documented here. However it seems that the tool according to this discussion was removed. And this bug report shows it was purged and is no longer distributed with QGIS.