texts: ArcSDE - Registering table with Geography spatial column fails with error -138 (SQL Server 2008)

texts:
I have a table of point features I am trying to register with ArcSDE (9.3.1 on SQL Server 2008) using the following command:

sdelayer -o register -l ThirdPartyWellsSurface,Shape -e np -C OBJECTID,SDE -t GEOGRAPHY -R 1047 -i sde:sqlserver:hou-sqldev10 -s hou-sqldev10 -D sde -u user -p pass

Where 1047 is a valid SRID (in the SDE_spatial_references table) for this data. It is failing with the error:

Error: The given coordinate references are incompatible (-138).
Error: Cannot Create Layer.

Looks like this used to be a bug back in 9.3 that got fixed. I checked and I do not have any null instances in the table. Googling and the ESRI sites didn't turn up much. Running a describe on a similar dataset with SRID of 1047 yields:

C:\WINDOWS\system32>sdelayer -o describe_long -l thirdpartywellssurfacetest,shap
e -i sde:sqlserver:hou-sqldev10 -s hou-sqldev10 -D sde -u user -p pass

ArcSDE 9.3.1  for SQL Server Build 1632 Thu Feb 26 12:05:37  2009
Layer    Administration Utility
-----------------------------------------------------
Layer Description ....: <None>
Database              : SDE
Table Owner ..........: SWNGIS
Table Name ...........: THIRDPARTYWELLSSURFACETEST
Spatial Column .......: SHAPE
Layer Id .............: 121155
SRID .................: 1047
Minimum Shape Id .....: 1
Offset ...............:
  falsex:       -400.000000
  falsey:       -400.000000
System Units .........: 1145324610.666666
Z Offset..............:          0.000000
Z Units ..............:          1.000000
Measure Offset .......: <None>
Measure Units ........: <None>
XY Cluster Tolerance .:          0.000000055879
Spatial Index ........:
  parameter:    SPIDX_GRID,GRID0=0.0178906151191854,FULL
  exist:        Yes
  array form:   0.0178906151191854,0,0
Layer Envelope .......:
  minx:      -124.52074,        miny:        24.45000
  maxx:       -72.80923,        maxy:        48.99887
Entities .............: npc
Layer Type ...........: SDE
Creation Date ........: 07/28/10 16:57:27
I/O Mode .............: NORMAL
Autolocking ..........: Enabled
Precision.............: High
User Privileges ......: SELECT, UPDATE, INSERT, DELETE
Coordinate System ....: GEOGCS["GCS_North_American_1927",DATUM["D_North_American
_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UN
IT["Degree",0.0174532925199433]]

And I know all of my points lie within that layer envelope:

Any ideas on what "The given coordinate references are incompatible (-138)" could be telling me here?


texts: Is an "Apply" method required to save programmatical changes to the QGIS project file?

texts:
I'm trying to change the nodata settings on a large number of color index raster layers in my project. I want to swap the original nodata with one specified by the user.

The Python code I've written works kind of like it should. The Nodata I specify is applied in the map, but the original Nodata value also remains (leaving me with two Nodata values.)

The change is also reflected in the "Layer properties -> Transparency" window. However, when I save the project and re-open it, the changes made by the script are lost. If I manually do the nodata changes in the GUI, they are saved to the project file as they should.

Am I missing some sort of "Apply" method in my script? Or have I misunderstood that dataProvider() is the correct class to do this?

from PyQt4.QtGui import *

#Display user prompt for nodata
qid = QInputDialog()
title = "Enter nodata value"
label = "Nodata value: "
mode = QLineEdit.Normal
default = ""

#Get input from user
nodata, ok = QInputDialog.getText(qid, title, label, mode, default)
#Cast user nodata value to int
nodata = int(nodata)

li = iface.legendInterface()
#Get all visible layers
visiblelayers = [layer for layer in li.layers() if li.isLayerVisible(layer)]

#Iterate over all visible layers
for layer in visiblelayers:

  #Get data provider
  provider = layer.dataProvider()

  #Do not use source nodata value
  provider.setUseSrcNoDataValue(1,False)
  #Use nodata value from user
  provider.setNoDataValue(1,nodata)

  print "Nodata set to",nodata,"on",layer.name()
  layer.setCacheImage(None)
  layer.triggerRepaint()


texts: How to update polygon longitude if negative?

texts:
I have a polygon, with longitude in [-180, 180], with 2 parts around the dateline.
I want to make a map centered on the date line.
My plane is to:
1./ update the longitude, so that if lon < 0 , lon=lon+360
2./ dissolve the two part to get a unique polygon.

My problem is now 1./
I'm trying:

ogr2ogr nodateline\new_polygons.shp org_polygon.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry, 360, 0) FROM org_polygon WHERE X(geometry)<0"

and the command returns the following error:

Warning 6: Normalized/laundered field name: 'ShiftCoords(geometry, 360, 0)' to 'ShiftCoord'

EDIT: I also tried with ST_X(geometry) instead of X(geometry): same error.

The error only occur with the "WHERE" clause.
a./ what does the error mean?
b./ what would be the correct syntax?


texts: How to update raster pixel values based on another raster with conditions in QGIS?

texts:
I have 2 land cover rasters. I want to update parts of the first raster with values from the second using QGIS.

The workflow should be like this:

1. Select values equal to 30 and 34 in the first raster.

2. Check if there is a value at the corresponding pixel in the second raster.

3. If there is a value in the second raster and it is different from 0,
   change the value in the first raster to the value of the second raster.

4. If not, do nothing.

How should I proceed to achieve such updating tasks? I think raster calculator is not supporting conditional statements, so I'm not sure how I can do that.

test files used are here.


texts: Change Lable (Name of countries) according to Attribute values ArcMap

texts: I'm creating a world map using ArcMap. I want to style the label in such a way that their name should be displayed as the attribute value. For example if I'm using population data then the China's font should be largest then other


texts: symbology changes colour when harvested from WMS to online map

texts:
I'm harvesting a web map service and from my legend the symbology is pink but the same feature on the online map is green?
I've had a look at the GetCapabilities file and opened the legend png and it's pink. Does anyone have an idea as to why the colour is changing once visualised on my map?

Thanks


texts: GeoJSON data and XYZ (meter) coordinates

texts:
I'm relatively new to mapping.  I have some data which I would like to save using the GeoJson format and display on Google Maps (Javascript API v3).  This works for my shapes which are in [Lng,Lat] format. I also have some features (Line String) which are in [X,Y] format, where the units of X,Y are meters. I have an origin point which is in [Lng,Lat] format.

I've read a little bit about using Tiles for X,Y coordinates, but if I understand correctly, my X,Y units are pixels, not meters.

Currently I have been saving my X,Y data twice; once as X,Y and once as Lng,Lat (I use google maps geometry tools to translate the X,Y points to Lng,Lat).  Is this my only option, or can I save the data in X,Y [meters] with an origin and have it display properly on the map?

Also, the points are 'relatively' close for geo data. The highest X,Y values would be around 3000m.

Edit: Adding (truncated) Sample Data

XY(Z) data with origin.  (data is LineString, origin is Point)

{
  "type": "FeatureCollection",
  "features": [
    {
      "geometry": {
        "type": "LineString",
        "coordinates": [
            [
              717,
              1246.3812,
              0
            ],
            [
              703.1146,
              1227.8316,
              0
            ],
            [
              688.9295,
              1209.5377,
              0
            ],
            [
              838.4507,
              1426.9744,
              0
            ]
        ]
      },
      "type": "Feature"
    },
    {
      "geometry": {
        "type": "Point",
        "coordinates": [
          -121.76410900803,
          36.575784480854
        ]
      },
      "type": "Feature"
    }
  ]
}

Same data as above converted to Lng,Lat GeoJson:

    {
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [
            -121.7560874662,
            36.58698064498
          ],
          [
            -121.75624282811,
            36.586814021396
          ],
          [
            -121.75640154205,
            36.586649694817
          ],
          [
            -121.7560874662,
            36.58698064498
          ]
        ]
      }
    }
  ]
}

Clarification, I need the XYZ data for processing server-side (and I want to save those XYZ coordinates for future use), but I want to display the data back on the Google Map, do I need to convert it back to Lng,Lat (2nd data set shown) or can I display the XY(Z) data w/ Lng,Lat origin directly by referencing some specific CRS?


texts: I am using QField on an Android tablet, I would like to know if I could use a web mapping service like open street map as a base map. Without having the actual file stored on the tablet?


texts: Using Field Calculator to compare two fields in att table and returning if equal or not in ArcGIS for Desktop?

texts:
Using 10.2.2 I am attempting to compare two fields in an attribute table using the Field Calculator.  I'd like to return a value of '0' if they are equal and a value of '1' if they are not equal.

There are 600,000 records to compare and the overwhelming majority will end up being equal.  So, in essence this is an exercise to identify the minority that are not equal.


texts: Unable to print anything or write anything to file in python after importing processing

texts:
I've been trying to get standalone QGIS python scripting to work, but for some reason after I import processing, I cannot print anything to the console, nor can I write anything to file. My end goal is to create a library of standalone processing scripts, so being able to output information to the console is important to me.

I slightly altered the code from Error: algorithm not found, since the solution given gave me the error of "no module named processing".

My code is as follows:

#!/usr/bin/python
print "1"
import qgis
print "2"
from qgis.core import *
print "3"
from qgis.utils import *
utils.uninstallErrorHook()
print "4"
app = qgis.core.QgsApplication([], True)
print "5"
import processing
print "6"

class DummyInterface(object):
    def __init__(self):
        self.destCrs = None
    def __getattr__(self, *args, **kwargs):
        def dummy(*args, **kwargs):
            return DummyInterface()
        return dummy
    def __iter__(self):
        return self
    def next(self):
        raise StopIteration
    def layers(self):
        # simulate iface.legendInterface().layers()
        return qgis.core.QgsMapLayerRegistry.instance().mapLayers().values()
iface = DummyInterface()
plugin = processing.classFactory(iface)

f = open("/home/alex/output.txt", "a")

a_list = processing.alglist()
for item in a_list:
    f.write(str(item))

f.write(len(a_list))

f.close()

When run from the terminal, this is what I get:

alex@blackbook-bear: ~
$ python test.py
1
2
3
4
5
alex@blackbook-bear: ~
$ cat /home/alex/output.txt
cat: /home/alex/output.txt: No such file or directory
alex@blackbook-bear: ~
$


texts: Using lat lon position to get data from polygon within shapefile using MapBox?

texts:
I want to use a dataset of historical counties and get a parish and county based on a lat lon position.

An extract of the data is here:

https://www.dropbox.com/s/ompf9nqjegcv1ts/par1851.csv?dl=0

The shape data looks like this:

So the idea is someone selects point on a map (Mapbox), this returns a lat, lon which I then want to use to dynamically return a 'county' and 'parish' based on the data in the csv.

I will need to use a system to determine if the lat, lon is within the area of a parish.

Does anyone has some pointers to get me started?

UPDATE:
Created example based on Bill's answer.


texts: List of coordinates for earth locations in Cesium?

texts:
All of the examples in Cesium show coordinates for the locations in the code. For example, for Wyoming they show the following. Where can we get these coordinates from? If I want to do Arkansas, is there a list of US states and their Cesium coordinates somewhere?

var wyoming = viewer.entities.add({
    name : 'Wyoming',
    polygon : {
        hierarchy : Cesium.Cartesian3.fromDegreesArray([
            -109.080842,45.002073,
            -105.91517,45.002073,
            -104.058488,44.996596,
            -104.053011,43.002989,
            -104.053011,41.003906,
            -105.728954,40.998429,
            -107.919731,41.003906,
            -109.04798,40.998429,
            -111.047063,40.998429,
            -111.047063,42.000709,
            -111.047063,44.476286,
            -111.05254,45.002073
        ]),
        material : Cesium.Color.RED.withAlpha(0.5),
        outline : true,
        outlineColor : Cesium.Color.BLACK
    }
});

viewer.zoomTo(wyoming);


texts: I'm having an issue making a figure for a poster with ArcMap.  I'm trying to make a grid that follows the map boundary in a Lambert Conical projection, and it's not coming out right.  The map boundary is a layer in my data frame, and I'm trying to clip the grid to that.  Everything except the bottom grid is working properly and I can't figure out how to fix it.  I've attached a screenshot for reference.


texts: Order of layers in PyQGIS layout function

texts:
I'm creating a standalone PyQGIS application with a function to print a layout with all layers on the existing canvas.
Here is the current code for printing all the layers on canvas onto a layout, and then exporting them into a pdf:
def createLayout(Main):
    layers = Main.canvas.layers()
    project = QgsProject.instance()
    manager = project.layoutManager()
    layoutName = 'Layout1'
    layouts_list = manager.printLayouts()
    for layout in layouts_list:
        if layout.name() == layoutName:
            manager.removeLayout(layout)
    layout = QgsPrintLayout(project)
    layout.initializeDefaults()
    layout.setName(layoutName)
    manager.addLayout(layout)

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

    ms = QgsMapSettings()
    ms.setLayers(layers)
    print(ms.layerIds())
    print(layers)
    rect = QgsRectangle(ms.fullExtent())
    rect.scale(1.0)
    ms.setExtent(Main.canvas.extent())
    map.setExtent(Main.canvas.extent())
    layout.addLayoutItem(map)

    map.attemptMove(QgsLayoutPoint(5, 20, QgsUnitTypes.LayoutMillimeters))
    map.attemptResize(QgsLayoutSize(180, 180, QgsUnitTypes.LayoutMillimeters))

    # This exports a Print Layout as an image
    layout = manager.layoutByName(layoutName)

    # This creates a QgsLayoutExporter object
    exporter = QgsLayoutExporter(layout)

    # This exports a pdf of the layout object
    exporter.exportToPdf('C:\\Users\\seteg\\Desktop\\test\\testPDF.pdf', QgsLayoutExporter.PdfExportSettings())

I was following this tutorial https://www.youtube.com/watch?v=h8tNuzHqK2I and it works for the most part. Problem is, in my canvas, there should be two layers, an shp layer and a tiff layer. The shp file should be on top of the tiff file.
This is the output of print(layers):
[<QgsMapLayer: '0-Boundary' (ogr)>, <QgsMapLayer: '1-1-Slope' (gdal)>]

As mentioned in the API for setLayers:(https://qgis.org/pyqgis/3.0/core/Map/QgsMapSettings.html#qgis.core.QgsMapSettings.setLayers)
"Set list of layers for map rendering. The layers must be registered in QgsProject. The layers are stored in the reverse order of how they are rendered (layer with index 0 will be on top)"
My Boundary Layer is index 0 while Slope layer is index 1 but this is my output:

The tiff layer is on top of the boundary layer.
What is the problem here? I have tried reversing the order of layers by replacing layers with layers[::-1] but to no avail.

EDIT:
running print(ms.layerIds())
give me
'0_Boundary_8deaf683_25ea_4413_ba85_9e932b10db32', '1_1_Slope_3b8789a9_8e2c_4d0f_a81f_ce80955bd4ff']

So I do not know why slope layer is on top of boundary.


texts: Is there a shortcut in GRASS for importing vector files?

texts: Is there a keyboard or command shortcut in GRASS for importing vector files? If you know link where I can find maybe whole list of shortcuts for GRASS, that would be awesome.


texts: How to create 3D web application using Open Source GIS software?

texts:
I'm using OpenLayers - 2.11,GeoServer 2.13, ExtJS 3.4, and Java for developing GIS based web applications.
How can we create GIS based 3d application like below attachment using opensource softwares in GIS.


texts: Using zonal stats on lidar data returns unusual values

texts:
I have a raster of tree heights from Lidar and a polygon layer where I have sketched the tree crowns by eye from an orthophoto.  See screenshot.  I want to get the maximum heights within each tree polygon so used zonal statistics as table, but the table returned a load of nonsense.  Pixel inspector on the raster layer shows sensible values (around 9 metres for small trees up to 23 metres for the bigger trees, which correlate well with field estimates) but none of the columns in the stats table (below) make sense i.e. the max value ought to equal the max values found with pixel inspector but they dont and the ranges are really big.
Any idea what could be the issue?


texts: Find feature by property from multiple layers

texts:
What is the best way to find a particular feature by a property (let's say id), where the feature could exist in one of many layers? I could theoretically step through all the layers, and step through all the features in those layers until I find it, but I was hoping there would be a better/faster way.

Basically, here is the way I'm doing it, but would rather a better way if one exists:

map.getLayers().forEach(function (lyr) {
  const layerSource = lyr.getSource();
  const features = layerSource.getFeatures();
    for (const feature of features) {
        if (feature.get("id")=== valueIAmLookingFor){
            //do something
        }
    }
});


texts: Labels misrendered in this Mapbox Studio Classic basemap

texts:
I created a base map with Mapbox Studio Classic, and I'm seeing a problem with how the street labels are rendering at close levels of zoom. For an example, I've highlighted (in yellow) some mistakes in the screenshot below. In some cases the label is just cut off in the middle; in others pieces of the label text are repeated inappropriately. (The two northwest street labels should read "Oakland" and "Maryland".)

I used the Light template, and the only change I made to it was to color the water blue. The mistakes appear in my original rendering within Mapbox Studio Classic, as well as in CartoDB, where I've loaded the tiles as a custom base map. Reloadings and zooming in and out do not fix the problem -- when I return to this zoom level, the mistakes are still there. And there are many throughout the map, not just the ones in the screenshot.

Is Mapbox Studio Classic still be supported? I'm not using the new Studio for this project because I need to be able to use the base map in Cartodb, which doesn't work with the new version.


texts: Label Position Changed In VectorTile Service

texts: I made a map of boundies in ArcGIS Pro, labeled the polygon geometry by it's name. Everything is okay in ArcGIS Pro, but When I publish it to the ArcGIS Portal as a VectorTile Layer some of the label positions change.


texts: Joining selected fields from table to shapefile using ArcGIS Pro?

texts:
I have a line shapefile with 100+ fields and a separate dbf file with 100+ fields.  I'd like to programatically (using arcpy and Python 3) eliminate all but 3 fields in both shapefile and the dbf, and then join on a common field.  The desired output would only contain 5 fields.

If I do a simple join:

joined = arcpy.AddJoin(path_to_shp, "id_field", path_to_dbf, "id_field")

I cannot delete the extra fields. And if I write this joined table to a file, the joined field names become unrecognizable due to length constraints.

I was looking at the 'field_info' parameter of MakeFeatureLayer, but I can't find any examples of how this would be implemented, and it seems like none of the tools I'd normally use for this take FieldMappings.


texts: Automatically georeference hundreds of png files based on polygon vector coordinates

texts:
I have 1000 polygons in my shapefile. I also have 1000 PNG files. Do you think it is possible (based on a common attribute), through a script or something, to assign each of those PNG files to a polygon? I am looking for an automated process which should "read" the polygon coordinates and based on that common attribute, it should somehow take the corresponding PNG, transform it in a tiff file and place it between the boundaries of that polygon? Basically, do the georeferencing for me ...

I use mostly open source software and Global Mapper but I have no clue on how to do this..

Many thanks for the support,


texts: Using ArcGIS, ArcInfo and fme, I would like to transform data which contains wkt geometries (points, lines and polygons) into SHP. Has anyone a workflow for doing such a transformation?


texts:
Where can I download the last default.style for osm2pgsql?

I get http://svn.openstreetmap.org/applications/utils/export/osm2pgsql/default.style - this directory does not exist when I try to download it.


texts: What data type is excel in Arcgis Add Script wizard?

texts: When creating tool from a python script in ArcGIS, data type has to be chosen in Add Script wizard and one of my parameters is Excel and another one specific sheet whereas there can be more sheets in excel. I don't know which one to choose for excel and sheet which is obtained from excel. Both excel and sheet are parameters in my script.


texts: How can I apply a style rule to multiple categories? (2.14 and above)

texts:
I think this used to be more straightforward in earlier QGIS versions, but with more powerful rule-based styling (great development) I am now a little confused as to how to apply a rule to multiple categories at once in the Style dialogue.

In the example below, I've styled my data using categories. I want to apply the $id = $atlasfeatureid rule to all categories because I'm working on an Atlas.

I'm guessing I'm missing something simple, but can't think what. It's not a problem to apply rule separately to a few categories but when working with 50+ it would be very time consuming.


texts: Error: Transforming SpatialPolygons to the CRS of the Raster

texts:
I want to get the mean of all raster cells in each part of a polygon but I am getting a warning message. My code is:
rain <- readOGR(dsn="/Users___.shp", layer="rain")
class(rain)
#SpatialPolygonsDataFrame

region <- readOGR(dsn="/Users___.shp", layer="region")
class(region)
#SpatialPolygonsDataFrame

s <- raster(ncol=90, nrow=90)
extent(s) <- extent(rain)
rain_raster <- rasterize(rain, s)
class(rain_raster)
#RasterLayer

projectRaster(rain_raster, crs = projection(region))
test1 <- extract(rain_raster, region, fun="mean", na.rm=TRUE)
#warning message: Transforming SpatialPolygons to the CRS of the Raster

How can I solve for this problem?


texts: Changing symbology at different scales in ArcMap?

texts:
I am using ArcMap 10 and I want to display my village parcel symbology differently at given scales. I have one idea to copy and paste same layer.

But my question is that in a single layer, is it possible to give different symbology as per the scale? e.g

Layer Name    symbol   scale
xyz            -----   100000
xyz            *****   200000


texts: How to download older version of arcgis?

texts: I need to re Install ArcGis after formatting one computer, but on version 10.1 (not 10.2.2 that is available on line). were can i found this version for download ?


texts: Dojo Promise/All Array Not Populating InfoTemplate in ArcGIS Javascript

texts:
Using dojo/promise/all I'm attempting to populate an infoTemplate with results from a queryTask which runs later in the script.

Our app returns an XY value from a geocoded address search and queries that value to see if it falls within a series of polygonal features. If it does it should pass a string for each of the three queries.

The query runs and the console doesn't log an error. However, the infoTemplate is blank:

I've tried to loop through the array results prior to instantiating the infoTemplate but I don't know that that's necessary. Also, I've found that I'm unable to manually set the title of the InfoTemplate using a string.

Code for adding graphic/calling array:

          function addPlaceGraphic(item, symbol) {
              map.graphics.clear();
              var place = {};
              var attributes, infoTemplate, pt, graphic;
              pt = item.feature.geometry;
              place.address = item.name;
              place.score = item.feature.attributes.score;

              var trashQuery = runPWQuery(pt, trashLayer, "trash");
              var yardQuery = runPWQuery(pt, yardWasteLayer, "yardwaste");
              var recycleQuery = runPWQuery(pt, recyclingLayer, "recycling");

              all([trashQuery, yardQuery, recycleQuery]).then(function (results) {

                  //receive query responses, update attributes and add graphics

                  attributes = { trashQuery: pt.trashLayer, recycleQuery: pt.recyclingLayer, yardQuery: pt.yardWasteLayer };

                  infoTemplate = new InfoTemplate();
                  infoTemplate.setTitle("Scheduling Information");
                  infoTemplate.setContent("Trash: ${trashQuery}
Recycling: ${recycleQuery}
Mixed Bulk: ${yardQuery}");
                  graphic = new Graphic(pt, symbol, attributes, infoTemplate);

                  //add to map
                  map.graphics.add(graphic);
                  map.centerAt(pt);
              });
          }

Code later in script for queryTask/returning promise:

function runPWQuery(in_geometry , in_fl , in_container_id) {
        var query = new esri.tasks.Query();
        query.returnGeometry = true;
        query.outFields = ["*"];
        query.geometry = in_geometry;

       var promise = in_fl.queryFeatures(query, function(myresponse, io) {
              var temp_val;
              var values = [];
              var tstr;
              for (var il = 0; il < myresponse.features.length; il++) {
                  if (myresponse.features[il].attributes["MONDAY"] == "Yes") {
                      temp_val = "Monday";
                  }
                  else if (myresponse.features[il].attributes["TUESDAY"] == "Yes") {
                      temp_val = "Tuesday";
                  }
                  else if (myresponse.features[il].attributes["WEDNESDAY"] == "Yes") {
                      temp_val = "Wednesday";
                  }
                  else if (myresponse.features[il].attributes["THURSDAY"] == "Yes") {
                      temp_val = "Thursday";
                  }
                  else if (myresponse.features[il].attributes["FRIDAY"] == "Yes") {
                      temp_val = "Friday";
                  }
                  else {
                      //temp_val = "Other";
                      temp_val = myresponse.features[il].attributes["DESCRIPT"];
                  }
              }

              var statCount = myresponse.features.length;
              if (statCount >= 1) {
                  $("#" + in_container_id).html(temp_val);
                                   }
              else {
                  $("#" + in_container_id).html("");
              }

              return temp_val;

          }, function (error) {
              console.log(dojo.toJson(error, true));
          });

          //return promise
          return promise;

    }

    function renderPWQuery() {

    }


texts: Rasterio creating composite/stack much larger than the individual bands

texts:
I am using rasterio to create a stack of eight bands and each individual band is about 1GB, but when it's stacking them the output (with a lzw compression) is a tif that's over 53GB! This seems huge to me. Anyone know if this is normal?

def stack_bands(list_of_rasterio_bands, temp_dir, outputname, output_ext):
    # this will stack bands in the order they are in the list, returns a rasterio object of what was created
    print("[INFO] creating {}".format(outputname))

    composite_path = os.path.join(temp_dir, outputname+output_ext)

    arr_stack = []
    for band in list_of_rasterio_bands:
        band_arr = band.read(1, masked=True)
        arr_stack.append(band_arr)

    out_meta = list_of_rasterio_bands[0].meta.copy()
    out_meta.update({"count": len(list_of_rasterio_bands),
                     "nodata": -10000,
                     "compress": 'lzw'})

    with rasterio.open(composite_path, "w", **out_meta, BIGTIFF="IF_SAFER") as dest:
        for band_nr, src in enumerate(arr_stack, start=1):
            dest.write(src, band_nr)
    print("[INFO] {} saved".format(outputname))

    return rasterio.open(composite_path)


texts: Copy multiple shapefiles that intersect a polygon to a new location

texts:
I am trying to select and copy shapefiles from one folder to a new location based on if they intersect a polygon feature class. I am getting stuck in somewhere in the for loop.
EDIT: The goal is to copy the whole file if the extent of the shapefile intersects the polygon. I am not wanting to copy a subset of any shape.
# Generate list of shapefiles in the workspace
# Copy shapefiles in a list that intersect a polygon to a new location

import arcpy
import os

# set environment settings
arcpy.env.workspace = 'E:\\workspace\\Contour_Shapefile'

# set local variables
out_workspace = 'C:\\localvariables\\outworkspace\\'
select_features = 'C:\\localvariables\\Project.gdb\\AreaOfInterest'

# use ListFeatureClasses to generate a list of shapefiles in the workspace shown above.
fc_list = arcpy.ListFeatureClasses()

#Execute CopyFeatures for each input shapefile that INTERSECT Area of interest.
i=1
for shapefile in fc_list:
    out_featureclass = os.path.join(out_workspace, os.path.splitext(shapefile)[0])
    arcpy.SelectLayerByLocation_management(shapefile, 'INTERSECT', select_features)
    arcpy.CopyFeatures_management(shapefile, out_featureclass)
    print(str(i) + ' copied ' + shapefile)
    i += 1
else:
    print( str(i) + ' MISSING')
    i += 1
inpfile.close()


texts: Map algebra with loop and SetNull not working

texts:
I'm trying to run a loop through a number of rasters while using SetNull but keep receiving an error message. I have worked on these rasters individually using SetNull and it it works, but would need to do this multiple times so a loop would be so much better.

I am new to Python.

Update - sorted it out! The script below shows the changes in comments

# basic mask example: SetNull according to criteria and output set to 1
# Description: Returns NoData to any input cell with a value greater than 120 in the output raster, and the remaining output cells will be set to 1.
# Requirements: Spatial Analyst extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = r'C:\Users\Documents\ArcGIS\pythontest\acc'
arcpy.env.scratchWorkspace = r'C:\Users\Documents\ArcGIS\pythontest\acc'

#Set local variables
rasters = arcpy.ListRasters()
#inRaster = rasters
inFalseRaster = 1
whereClause = "value > 120"
#output folder is now different bc don't want to overwrite
outfolder = r'C:\Users\Documents\ArcGIS\pythontest\test'

#Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Loop through list of rasters
for raster in rasters:
    #try/except has been dropped
    # Execute SetNull outSetNull = SetNull(inRaster, inFalseRaster, whereClause)
    #outSetNull = arcpy.gp.SetNull_sa(inRaster, inFalseRaster, "value > 120")
    # updated with correct variable from the for loop
    # rephrased the SetNull function. arcpy.gp.SetNull_sa is not correct
    outSetNull = SetNull(raster, inFalseRaster, "value > 120")
    #print the output before attempting to save
    print(outSetNull)
    #import os module
    import os
    #Set output name for each output to be the same as the input
    # updated with correct variable from the for loop, attempting os module
    outSetNull = os.path.join(outfolder, raster)
    outSetNull.save(outras)

Old error message - refers to the line where    outSetNull = arcpy.gp.SetNull_sa(raster, inFalseRaster, "value > 120") :

Runtime error  Traceback (most recent call last):   File "", line 31, in    File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing_base.py", line 510, in      return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 000877: Output raster: The grid or stack base name in C:\Users\Documents\ArcGIS\pythontest\acc\value > 120 has space.


texts: How do you symbolise a polygon to look like an AutoCAD revision cloud?

texts:
I work for an engineering firm and the standard for annotating areas which require review or areas of interest is to add AutoCAD revision clouds. There is a dedicated tool for this in AutoCAD but I want to create them in QGIS.

Drawing them freehand will take too long so is there a way of symbolising a simple polygon to look like a revision cloud or alternatively is there a way of drawing a polygon which will automatically take on the shape of a revision cloud?

AutoCAD revision cloud examples:


texts: QGIS Raster file not matching with vector data when I match their Coordinate referencing systems

texts:
I'm still quite new to GIS. I am trying to sync the CRS of some vector data I retrieved in the field with some raster Digital Elevation Models I have downloaded. The vector data is a gpx. file referenced to WGS 84. The raster images do not appear to have an inherent referencing specification.

When I assign both to WGS84 it doesn't work. I have tried reprojecting the raster data.

I have tried some different things, but my vector and raster data are appearing in entirely different areas of the map field. When trying to assign both to WGS 84, the coordinates of my vector data is: WGS: 151.93694,0.00000. The coordinates of my raster data are about:86.427034,0.000000
I have also tried to set both the vector and raster files to GDA94 referencing and when reprojecting the vector data I ran an error:

Thanks for any help you can offer and pls let me know if you would like more specific info.
UPDATE: this is the transrmfoation prompt that comes up when I add the raster (GDA94) files to the project after having added the vector (WGS86) data with the CRS set to 'default to layer'.


texts: Using attribute-based line to style lines with different Color for WMS in GeoServer?

texts:
I am trying to use the styling of WMS using the so called SLD. I never used it before and I was wondering if I can do styling of a layer which has lines based on property.
I found here that something like this is possible. They call it (attribute-based line).
What I don't understand is where to define these attributes.
My line data are stored in a postgis and loaded in a store in geoserver. My postgis table has 3 columns (id, geometry, typeOfRoad).
My question is: Can I use the typeOfRoad attribute to color different the lines different?
Is there an example?


texts: Using GeoPandas, Shapely and Fiona to "select" state from shapefile, save it and display it?

texts:
I have the usual US census states shapefile for the 48 contiguous states. I am trying to select a state of interest from that shapefile then save that shapefile and display. I want to use GeoPandas, Shapely, Fiona to do this.

Any ideas on how this could be done?


texts: Rasterizing vector data for landscape permeability analysis

texts:
So I vectorized an area with 12 different categories/features. My intention is to use this in the future for calculating landscape permeability.

I would like to know which steps I need to take right now.
I presume it would be better to transform this into a raster, right? So I can give each pixel a value (in this case, 12 different values) and after that, continue with r.walk and so on.

If so, do I need to merge all these 12 files into one? Can I transform this into a raster while still keeping each of the different categories? (including the line vectors?)


texts: Intersection of two vector tables data in postgresql

texts:
I have satellite raster data that is converted into pixels using

CREATE TABLE gtru_tile_pix AS SELECT
    (ST_PixelAsPolygons(rast)).val,
    (ST_PixelAsPolygons(rast)).geom
FROM gtru_tile_raster;

and I have vector data (polygons) in other layer to get sample pixels from the above pixels.

Please tell me
How to get the intersection of two tables?


texts: Adding values of common field to field during union?

texts:
I have 3 overlay features; A,B,C and each has a field called FIELD_ID

If I do a union, I will get 3 new fields in the output: FIELD_ID_1, FIELD_ID_2, FIELD_ID_3

But how do I get just one summary field 'FIELD_SUMMARY' that contains the values of FIELD_ID_1, FIELD_ID_2, FIELD_ID_3 for each polygon?

How do I list the values of [FIELD_ID] for features A,B and C  that occur from the union but in a single field, not multiple ones?

I'm aware that I can do a field calculation after the union that will concatenate values from FIELD_ID_1, FIELD_ID_2, FIELD_ID_3 etc, but that seems a bit long winded


texts: Using dynamic text to display layer credits

texts:
I'm looking for a way to update the 'layer credits' into a dynamic text box. I've found ways of using the map credits or data frame credits here but can't seem a way to get the information derived from a layer.

I just seem to produce an empty string every time. Has anyone done this or found a work around? Seems like an extremely useful way to update map credits with little to no resources.


texts: Add a list of dates as an image properties in Google Earth Engine?

texts:
I've merged two image collections after calculating median values by month in each collection from 2000 to 2018. However I need to carry out an operation that needs the "sytem:time.start" field that was lost after calculating the monthly means. Therefore, I would like to know how to add a property within each image of the collection from a list representing month and year? The only property I have now in my collection is "system:index".

The code can be accessed by this link.
The collection I want to associate with "system:time.start" information is on line 164, named "joined_collections"


texts: converting vector points to lines using QGIS 2.8

texts:
I have about 350k of vector points that represent parking spaces, I need these points to be converted to lines and connected to the nearest point available.

I've tried the point2one plugin, but connects parking points that aren't relevant to each other, and I've tried doing a nearest neighbor analysis in GRASS, but the data is to large and just crashes the program.

any suggestions?


texts:
I'm trying to make a survey of all regions in Sweden for an article that me and my colleague are working on. When I'm trying to import the shapefiles into CartoDB it says " Unfortunately this file extension is not allowed ". Im guessing the reason is that i got the free version.

Is it possible to work this out in some way?


texts: Problem when inserting a QLR file in QGIS

texts:
I work with a polyline layer:

and the layer located in "G:\desktop\qgis\lyr".

After i save the layer in qlr file and saved it in "G:\desktop\qgis\lyr", i try to load the qlr file with "add from layer definition file" function,to a map that located in "G:\desktop\Project".

The result is that the qlr file doesn't appear on the map and there no properties for that layer. There no symbology either.


texts: Arcpy Update Cursor-Overwrite values in field based on CSV

texts:
I have a CSV file called Modifications. The header has two entries, Rental Object and Rental Object Name. The Rental Object is a string that has matching values in a FGDB feature class where the column is called Property_Ref. This attribute table also contains a field called Descriptio. I need to update this Descriptio with the values from Rental Object Name in the CSV. The Rental Object in the CSV and the Property_Ref match. I'm not sure how to complete this. Below is as far as I've progressed.

# Update Records to be Archived:
import csv, arcpy
fc = "U:/Data/Lease_Plans.gdb/LeasePlans_SAP_NOs_Descs"
Modifications= "U:/Data/Rental_Object_Modifications.csv"

#Read the CSV File
with open(Modifications, 'r') as f:
    reader = csv.reader(f)
    mydict = dict((rows[0],rows[1]) for rows in reader)

fc_fields = ('Property_Ref','Descriptio')

with arcpy.da.UpdateCursor(fc, fc_fields) as cursor:
    for row in cursor:
        #Not sure how to structure this.

del cursor


texts: Is it possible to create one single NDVI image out of multiple NDVI images?

texts:
I am working with Landsat 8 images to create seasonal NDVI images in QGIS.

For all NDVI images I masked the clouds/cloud shadows so there are 4 different NDVI images with different cloud/cloud shadow masks .

How could I create one single NDVI image with mean NDVI values out of the 4 different NDVI images?


texts: A simple label for an entire MapServer mapfile?

texts:
I want to incorporate a simple text label into a mapfile-based map, such that the file's resulting images always display a specific text string (ideally in white over a black background) along the top of the image.

The image's pixel dimensions are known, but the map's extent will depend on outside factors; I'd rather not orient the label to spatial coordinates, but if I have to I will.

Despite how simple this seems to me, I can't seem to find any way to do it, including with legends (which don't seem to be super-imposed on the map anyway, but somehow work differently).


texts: QGIS gets stuck when running script to automatically refresh layer

texts:
I have a QGIS layer that represents a set of links and nodes to show the streets of a city and the traffic speeds on each street. I've been trying to write a script so that if I change the speeds in the CSV file (links3.csv) the speeds will automatically update in QGIS. However when I run the code I currently have, then change the file, the entire program freezes and I have to force quit it.

import qgis.utils
from PyQt4.QtCore import QFileSystemWatcher

myLayer = qgis.utils.iface.activeLayer()
def refreshLayer():
    myLayer.setCacheImage( None )
    myLayer.triggerRepaint()

watcher = QFileSystemWatcher()
watcher.addPath( '/Users/username/Traffic Speeds/links3.csv' )
watcher.fileChanged.connect( refreshLayer )

What am I doing wrong? I got the code to automatically refresh the layer from here


texts: Display GeoServer's UTFGrid output in OpenLayers 3

texts:
GeoServer added UTFGrid output for vector features starting from the 2.9-beta version. It is implemented as a WMS service output format.

After installing GeoServer 2.9-beta, adding a datastore with a PostGIS connection and publishing a layer that contains vector features, I can click the WMS UTFGrid entry from the available formats list on the Layer Preview page in the admin UI and indeed get a UTFGrid representation of my data displayed in the browser.

How can I now integrate this as a UTFGrid layer into OpenLayers 3? The only resource that I could find on this topic is this official example. However, it uses ol.source.TileUTFGrid as the layer's source and passes a URL which ends in .json. GeoServer, however, publishes a WMS URL which expects further parameters such as request, bbox etc.

Simply passing the URL generated by the Layer Preview page gives Uncaught SyntaxError: Unexpected token : in the browser's console.

What would be the correct way to include GeoServer's UTFGrid output as a layer?


texts:
I have followed the instruction on the web page

https://faunaliagis.wordpress.com/2013/08/15/adding-ecw-support-in-qgis/

to add ECW support to QGIS but it seems to have no effect on QGIS server. I cannot see the raster layers in a QGIS project when I query the QGIS server capabilities. I have also noticed that the raster layers are not listed when using some other raster formats, but they are listed if I use TIFF, or IMG format.


texts:
I have a layer of polygons, cadastrals, that i'd like to have the length of the edge of:

I thought about doing this:

Dissolve the polygons
Create one large dissolved polygon
Making a; From polygon to Line

But from here, i don't know what to do next.

Cut it up, and measure the length of the line, then join attributes so that the polygons get the data. But how?


texts:
I am newbie in Geomajas and I tried to create new Geomajas application using the documentation in http://files.geomajas.org/maven/trunk/geomajas/docbook-gettingstarted/html/master.html.

I don't know how to display my country's map (Mongolia) in my created Geomajas application. Can you tell me how to do that or do you have any step by step guide/tutorial on this issue?

Please help me


texts:
At 3.12.3, I need to group the two items highlighted in this screenshot:

According to the QGIS documentation: "To create a group of items, select more than one and press the Group button..." [Emphasis added]
But HOW do I select more than one item???
Here's what I've tried:
In my example above, I select the first item by clicking on either one of the two, which highlights it.  I then ctrl-click the second item.  But that just switches the selection and highlights the second item; it is not added to the selection.  Same thing when I try shift-click.
Thanks in advance!


texts: Load support data set from NetCDF file (Tropomi Sentinel 5P) with GDAL.NET

texts:
Been trying to load two data sets from a Tropomi NetCDF file (N2O measurements). I can manage to load some data sets, which are right under PRODUCT, but not others which are under PRODUCT/SUPPORT_DATA/GEOLOCATIONS/.

I've been trying to load them the same way as I do the other datasets.

using (Dataset ds = Gdal.Open(filePath, Access.GA_ReadOnly))
{
  list = ds.GetMetadata("SUBDATASETS").Cast<string>().ToList();
}

The first thing that draws attention is that the PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds and PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude bounds do not appear in the SUBDATASETS.

If I try to address it manually, I do get a dataset back, but it's invalid as any subsequent read from it fails:

So, this fails:

var latBounds = Gdal.Open("HDF5:D:\\Downloads\\tropomi_no2_20181227\\S5P_OFFL_L2__NO2____20181227T114003_20181227T132133_06244_01_010202_20190102T140322.nc://PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds", Access.GA_ReadOnly);
float[] latBoundsBuffer = new float[latBounds.RasterXSize * 1];
var latBoundsRaster = latBounds.GetRasterBand(1);
var scanline = latBoundsRaster.ReadRaster(0, 0, latBounds.RasterXSize, 1, latBoundsBuffer, latBounds.RasterXSize, 1, 0, 0);

Any idea how these auxiliary datasets can be loaded? If useful, the product data sets are described here.

Also, Panoply shows and load the datasets just fine:


texts: Raster masking with vector error in QGIS 3.6

texts:
I have a global glc-2000 raster map. I would like to clip it to study specific area. I used boundries of continents shapefile as mask. Both of them have same projection.  Since I am beginner I have followed steps of this tutorial using QGIS 3.6:
http://www.qgistutorials.com/en/docs/raster_mosaicing_and_clipping.html

I  hoped that I get Europe output raster layer but the result was a black hole. What went wrong?


texts: Passing where clause as variable to arcpy.SearchCursor() gives ERROR 999999?

texts:
I'm having trouble passing my where clause for arcpy.SearchCursor as a variable.

The following code successfully generates a query for me:

PolyList = "C:\\Users\\korora912\\Desktop\\polygons.shp"
cursor = arcpy.SearchCursor(PolyList,"""date > date '2016-01-01 00:00:00' AND date < date '2016-01-04 00:00:00'""")

However when I try to pass the dates as variables, like so:

StartDate = "2016/01/01"
date_1 = datetime.datetime.strptime(StartDate, "%Y/%m/%d")
end_date = date_1 + datetime.timedelta(days=3)
where_clause = "\"\"\"date > date '%s' AND date < date '%s'\"\"\"" %(date_1,end_date)
cursor = arcpy.SearchCursor(PolyList,where_clause)

I get: "RuntimeError: ERROR 999999: Error executing function.
An invalid SQL statement was used."

This, despite the fact that

print where_clause

returns the exact string that succeeds above:

"""date > date '2016-01-01 00:00:00' AND date < date '2016-01-04 00:00:00'"""

I know I must be doing something wrong, but with my limited Python and SQL experience, I'm at a loss. Can any of you with a more robust background help?


texts: How is scale for print composer calculated?

texts:
Not really an issue, more of a theoretical/technical question.

I'm doing regional maps with over 1500 km difference (~10 degrees of latitude) between lower/upper border  of the map, and I noticed that scale in print composer is fixed and does not change if you move scalebar within composer. So how is it calculated? Is it taking an average scale of the referenced map? A 100 km section(length) of the scalebar on the south would only equal ~80(?) km if moved to the north.


texts:
I've got a Python script for ArcGIS that I'm working on, and I'd like to have the ability to have the script quit if it doesn't have the necessary data available.  I tried a straight up sys.exit() but that would give an exception in ArcMap that I'd like to avoid.  I found This thread that suggests using a try block, so I made this function:

def quit_script(message):
log_msg(message) # already defined; writes a message to a file
if log_loc:
    output.close() # close the file used with log_msg()
try:
    sys.exit()
except SystemExit:
    pass

Unfortunately, that didn't work either.  Well, it doesn't make that error on ArcMap anymore, but it also doesn't, well, quit.  Right now, I have the bulk of my code in an if/else statement, but that's ugly.  Anybody have any other suggestions?

Thanks!
Brian


texts: Spatial join with concatenation of results

texts:
I've two polygon layers which I would make a spatial join.
The objective is to transport a specific code from layer A (Municipality) to layer B (Issues) using like a concatenate function, in this way I'll have in the layer B a field with the codes provided from the layer A. For instance if I've two polygon of A intersecting the same polygon on B I'll have a code like "A.B".
The problem is that ArcGIS with "spatial join" make this work but don't separate the value eg:AB instead A.B.

Some suggestions?


texts:
I'm trying to work on a school project, but I only need streams up stream from a certain point. As you can see in the picture, I've run the stream trace tool upstream from a point i've created.

I'm trying to find the total length of the red streams. However, you cannot select attributes from a stream trace. How can I accomplish this task?


texts: How can I export each polygon-vertice's coordinate (every point coordinate of every polygon) of a polygon Geodatabase-feature class to an ascii text file (csv or txt)? I don't want to convert it to a point FC first. I cannot use ET Geowizards! I don't have FME (or is the little FME encluded in ARC able to do this? How?) I don't want to export to XML format - except if someone can tell me how to convert XML to a simple text file.


texts: How to publish a postgis layer to geoserver?

texts:
Now my geoserver is on the amazon instance and the postgresql is on the other Amazon RDS. I have try the way provided by official document, but it always show the errors as following:
  "Caused by: org.apache.commons.dbcp.SQLNestedException: Cannot create
    PoolableConnectionFactory ..."


texts: Add-in installed, but its elements not appearing in ArcMap

texts:
I'm creating an Add-in for ArcMap 10.1.

When I export file testproject.esriaddin and install (whether double-click or from Customize Mode > Add From File dialog), message of success appears. But then, I can't find any of my elements in Customize Mode dialog (no Java Tools or Java Commands category, no toolbar in toolbars...).

I found this topic, which seemed to be about the same problem. There, the solution was to select JRE: Java SE 1.6 instead of 1.7. I set this in my project, but the Add-in is still not working.

What could be wrong?

Here is my java class for the button:

package cz.upol.geoinformatics.veselka;

import java.io.IOException;
import javax.swing.JOptionPane;

import com.esri.arcgis.addins.desktop.Button;
//import com.esri.arcgis.arcmapui.IMxDocument;
import com.esri.arcgis.framework.IApplication;
import com.esri.arcgis.interop.AutomationException;

public class Button1 extends Button {
    private IApplication app;
    //private IMxDocument projectDocument;

    /**
     * Called when the button is clicked.
     *
     * @exception java.io.IOException if there are interop problems.
     * @exception com.esri.arcgis.interop.AutomationException if the component throws an ArcObjects exception.
     */
    @Override
    public void onClick() throws IOException, AutomationException {
        JOptionPane.showMessageDialog(null, "Hello, World!");
    }

    @Override
    public void init(IApplication app){
        this.app = app;
        //this.projectDocument = (IMxDocument)app.getDocument();
    }
}

Here is how it looks after deploy to ArcMap:

I use ArcGIS + SDK 10.1 on Windows 7 32bit and Eclipse 4.3
Did I forgot to mention some details?

Edit: Here is my another screenshot: (You can see, that category Java Commnads is not listed in Categories)


texts: Seeking shapefiles of business locations

texts: Does anyone know where to find shapefiles of businesses such as grocery stores or gas stations?


texts: Adequate CRS when creating a grid over a whole continent in R

texts:
I'm trying to generate a grid of 6.72 km by 6.72 km points over the whole continent of Africa, but I'm having a lot of trouble in doing so.
I'm mainly running into 2 problems:
My first issue arises when I try to generate a grid on a single country. I downloaded the country level map of the world from GADM, kept only countries in Africa, and then I used the "Simplify Tool" In QGIS to make the polygons more simple so the code runs in a reasonable time.
Here is some example code for Angola based heavily on this answer on Stack Overflow, where I work with the original non-simplified GADM map, so I also changed the grid size to 200km so the code runs in a reasonable time. I had to change the CRS to EPSG:26910 for the code to run.
library(sf)
library(tidyverse)

#importing map of Angola
angola <- readRDS(gzcon(url("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_AGO_0_sf.rds")))
#checking shape
plot(angola$geometry)
#checking CRS
st_crs(angola)
#changing CRS
angola <- st_transform(angola, 26910)
st_crs(angola)

#building grid
grid <- angola %>%
  st_make_grid(cellsize =200000) %>%
  st_intersection(angola) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

#checking grid
plot(grid$geometry)

The problem here with the CRS is that the map with the grid is flipped. I tried changing back the CRS of the grid to the original after doing the calculations, but the dimensions of the pixels are not the desired when going back to the original CRS.
What CRS is the most adequate for this task? I'm very new to GIS, but I need a CRS that works properly in all Africa and that allows me to define a grid in meters and not in degrees.
My second issue is that I can't generate a grid over the whole continent of Africa. I downloaded the country level map of the world from GADM, kept only countries in Africa, and then I used the "Dissolve Tool" in QGIS to combine the shapes of all the countries and get the shape of the continent.
I can't provide code for people to replicate what I did in QGIS, so I attach this Dropbox link with the continent shapefile, and an example again with a grid with 200km, and that obviously suffer from the CRS problems mentioned in the previous point.
africa_comb <-st_read("africa_comb.shp")

#changing CRS
africa_comb <- st_transform(africa_comb, 26910)

#making the grid
grid_comb <- africa_comb %>%
  st_make_grid(cellsize =200000) %>%
  st_intersection(africa_comb) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

Here I get the following error:
Error in CPL_geos_op2(op, x, y) :
  Evaluation error: TopologyException: Input geom 1 is invalid: Self-intersection at or near point...

I have tried before doing the same without dissolving the countries, but I run into the same error. Both maps look like this after changing the CRS, so I'm suspecting the CRS might be causing the issue here:


texts: How do I stop the map at just one world/globe?

texts:
I am working on a project to put google maps layer and adding our data layers on top of it. using openlayers and google maps layer with out layer of top of it

now the google layer has max resolution of 156543.0339 and 22 zoom levels, I can limit the zoom levels and max resolution etc but how do I stop the map at just one world/globe

i.e. I don't want it to scrolling off to another wrapped around world!

Please see this fiddle for example,  http://jsfiddle.net/eaJHn/5

Anyone !!! ?


texts: How to rewrite a csv for QGIS to avoid "Invalid record format"?

texts:
So,
I have a .csv file.

Here an example how the a row from the original .csv looks like:

"AREA",44.37,-00.83,"5","1","8","E","E",,,"KKI","KKO"

Some lines have missing latitudes and longitudes, so I want to get rid of them.
Should be easy. Thats my snippet (Python 2.7) that does it:

import csv

n = 0
m = 0

nf = open("newfilefile.csv", "w")

with open("oldfile.csv", "r") as csvFile:
    originalData = csv.reader(csvFile, delimiter=",")
    for lines in originalData:

        if len(lines[9]) < 1:
            n += 1
            print 'Missing Latitude', n

        elif len(lines[10]) < 1:
            m += 1
            print 'Missing Longitude', m

        else:
             nf.write(str(lines[0]).strip() + "," + str(lines[1]).strip() + "," + str(lines[2]).strip()  + "," + str(lines[3]).strip() + "," + str(lines[4]).strip() + "," + str(lines[5]).strip() ...  + "\n" )

nf.close()

Thats how the same row would look like in my new .csv:

AREA,44.37,-00.83,5,1,8,E,E,,,KKI,KKO

Beforehand I could open the original csv file. QGis just said that there are missing geometry information. Now I get a bunch of errors, namely:
Invalid record format at line 324
and so forth... I compared the lines where it marks the error and the only difference is that in the oldfile, everything is wrapped in " while in the file my script creates that is not the case.

Why does QGis need " around every from a column (if that is the mistake) and how would I add that to my snippet?


texts:
When I start qgis I get a python exception:

Couldn't load plugin fTools due an error when calling its classFactory() method
...
309, in _import
   mod = _builtin_import(name, globals, locals, fromlist, level)
 File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/random.py", line    47, in
  from os import urandom as _urandom
ImportError: cannot import name urandom

From ImportError: cannot import name “urandom” [closed] it seems I would need to replace the python interpreter used by qgis but can't find a way to do that.

Any ideas on how to resolve / work around this?


texts: Use delimited text Vector data provider in QGIS 3.0 Processing Algorithm

texts:
I am writing a Processing Algorithm based on this template. Can I populate the QgsProcessingParameterFeatureSink from a delimitedtext provider in a similar fashion to the iface.addVectorLayer() (i.e. the following code snippet) ?

uri = ('file:///{}?'
       'type=regexp'
       '&delimiter={}'
       '&skipLines={}'
       '&useHeader=No'
       '&trimFields=Yes'
       '&xField=field_11'
       '&yField=field_12'
       '&crs={}'
       '&spatialIndex=Yes'
       '&subsetIndex=no'
       '&watchFile=Yes').format(myFile,
                                myRe,
                                countSkips(myFile),
                                myCrs)

layer = iface.addVectorLayer(uri, 'MY LAYER',  'delimitedtext')


texts: Group raster data by values in a second raster

texts:
I am looking for a PostGIS equivalent to the ArcGIS Zonal Statistics by Table function.

I have two rasters. I would like to group the data in the second raster by the values in the first raster, and report the summary statistics (count, sum, etc.) as a table. There are questions for grouping raster data by polygons, but I don't see anything relevant to grouping by rasters.

The desired result is produced by this query, but I am working with a global 30-arc second resolution raster, and would like to avoid dumping all the values. Note that the two rasters are aligned and the rids match for intersecting tiles. (I assume that unnest preserves the order of arrays - is this correct?)

SELECT id, sum(val), count(val) from (
   SELECT unnest(ST_DumpValues(r1.rast,1)) id,
          unnest(ST_DumpValues(r2.rast,1)) val
   FROM rast1 r1, rast2 r2 WHERE r1.rid=r2.rid) t1
GROUP BY id;


texts:
I've the following table in PostGIS

-- Table: public.elevation_data

-- DROP TABLE public.elevation_data;

CREATE TABLE public.elevation_data
(
  rid integer NOT NULL DEFAULT nextval('elevation_data_rid_seq'::regclass),
  rast raster,
  CONSTRAINT elevation_data_pkey PRIMARY KEY (rid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE public.elevation_data
  OWNER TO username;

-- Index: public.elevation_data_st_convexhull_idx

-- DROP INDEX public.elevation_data_st_convexhull_idx;

CREATE INDEX elevation_data_st_convexhull_idx
  ON public.elevation_data
  USING gist
  (st_convexhull(rast));

created with the following command

raster2pgsql -I -s 4326 -t auto DataSet\Elevation\*.dt2 public.elevation_data | psql -U username -d databasename

Now I want to export this table to a .dt2 file, so I can exchange it. How can I export the table in the raster format?

EDIT:

I've looking documentation and I've tried following query:

COPY(SELECT ST_AsGDALRaster(rast, 'DTED') FROM public.elevation_data) TO 'd:\dteddata.dt2' (FORMAT binary);

It creates a very big file (over 7Gb) when the table was created by importing two files of about 50Mb in total. When I try to open this file, for example in QGIS, it says me that the file is not supported. The same if I try to save it as GeoTiff with the command:

COPY(SELECT ST_AsGDALRaster(rast, 'GTiff') FROM public.elevation_data) TO 'd:\dteddata.tiff' (FORMAT binary);

even if in this case the file has a size of about 51Mb, that seems more correct to me.

I've followed instruction from this page for retrieving DTED data, and joining them with this other one, because ST_AsGDALRaster returns a bytea, and the second links explains how to save it to a file, but I'm doing it the wrong way of course...


texts: House address is not listed on Google Maps

texts:
When I search for my address on Google Maps it cannot be found which means I can't order takeaway etc.

Is there a way of adding my address to Google Maps or even add a street address?

I want the address to be publicly available so others can search it.


texts: What are Valid feature class and table names in ArcGIS Desktop?

texts: What are the rules for naming valid feature classes and tables in ArcGIS Desktop, and is there any supporting document which can explain all the dataset naming limitations?


texts: Openlayers 5 change kml style for icon only

texts:
I'm importing KML into my map where it has some point placemarks with images that are WAY too big.  What is the best way to import KML in the map, and just change 1 aspect of the style ( https://openlayers.org/en/latest/apidoc/module-ol_style_Icon-Icon.html , change the 'size' to a different value), but keep everything else?  For example, i'm adding the KML to the map in a very simple way:

        const kmlSource = new VectorSource({
          url: kmlURL,
          format: new KML()
        });
        const kmlVector = new VectorLayer({
          source: kmlSource
        });
        map.AddLayer(kmlVector);


texts:
As par the USGS page explaining the 8 bit quality image in following link:

http://landsat.usgs.gov/LandsatLookImages.php

Water is represented by bit 3 [Translated to color range (8 -15)] and Snow/ Ice is represented by bit 5 [Translated to color range (32 - 63)].

However, for my study location, the QA image shows water bodies as Snow/ Ice (Teal colored reaches) (color 32 -63) instead of (color 8 - 15) (Royal blue color) as par USGS reference. Any explanation?


texts: PyQGIS/PlugIn: How do I pre-check if processing function is executable before I get a QGIS critical error?

texts:
I built a PlugIn with the 'native:shortestpathpointtopoint' tool to find the shortest path between two points on a multipolyline. When I run the Plugin I want to get a message with information about wether the params are valid and the processing is actualliy feasable!

I can check the params for correct datatype, but how can i pre-check if the 'native:shortestpathpointtopoint' will actually work or not without getting a critical error. It wont work, when the two points are not connected by the polyline. I don't want a cryptic QGIS error but rather a: 'Pls check your points' messagebox.

I tried stuff like this but those are already pretty desperate tries:

def getPath(self):

    print('Pathfinder START')

    params = {
        'INPUT': lyr,
        'STRATEGY': 0,
        'TOLERANCE': 2,
        'START_POINT': self.Startpunkt,
        'END_POINT': self.Endpunkt,
        'OUTPUT': os.path.dirname(__file__) + '\\shortestPath\\shortestPath.shp'
    }
    if 'native:shortestpathpointtopoint'.canExecute():

        processing.run('native:shortestpathpointtopoint', params)

    else:
        print('Check your points')

    print('Pathfinder END')


texts: How to disable all mouse events in OpenLayers map div?

texts:
version: OpenLayers-2.13.1

brower : IE8 (I Hit it!!! But I need to work on it)

I create a div with id 'map' to display OpenLayers map. Sometimes, I need disable map div catch any mouse events, let me enable "drag" items into map (not really drag and drop event). But I don't know how to disable the map receives mouse events, so, I dynamically create a new div into the HTML body to cover above the map, make the new div receive all mouse events and without pass these events to the map, by this way, Google Chrome works well, but in IE8 is not!!!! So..... sad!!

In IE8, the created new div only receive mousemove event at itself's border, when mouse inner blank area of the new div, the mousemove events be received by the map div, upper level's new div receive nothing. In Chrome, when mouse move onto the new div, no matter on the border or inner blank area, all mouse event pass to upper level's new div first, the map div would not receive!

Could anybody told me how to make OpenLayers map div avoid receive mouse events and how to restore receive mouse events ability? If use new div cover upper the map, how to let upper div receive all mouse events in IE?


texts: Group PostGIS table-rows where geometries are close

texts:
In OSM there are often e.g. lines (paths) with the same name but different OSM-IDs. Right next to each other, so basically the same thing, but just different objects with no real relation to each other.
So in my table it would look like the following
osm_id    name      highway    geometry
---------------------------------------------
1         trail_1   path       <some-line>
2         trail_1   path       <some-line>
3         trail_2   path       <some-line>

I looked at ST_ClusterDBSCAN and that looks promising?!? I am just too bad with SQL to fully understand it.
So I'd like to group my lines by name, cluster them, look at each cluster -> take the metadata (osm_id, name, highway,... ) of one of them and also the combined geometry (ST_LineMerge or ST_Union?!?) of all of them -> Save as a new entry in a separate table.
So from my example in a separate table I only want to have TWO entries left:
osm_id    name      highway    geometry
---------------------------------------------
1         trail_1   path       <some-merged-line>
3         trail_2   path       <some-line>

I am not really sure yet how to tackle the whole thing. Can I do that solely in PostGIS?


texts: See WFS or WMS corresponding full http request in QGIS when a layer was added from any of these services

texts:
I have added both a new WMS and a new WFS server into a QGIS project, by specifying a generic GetCapabilities URL in the appropriate sections in the browser panel:

After that, I am able to add some of the available layers from these WMS and WFS servers.
Now, how can I find the full corresponding underlying http request (that was probably build under the hood by QGIS), so that I can copy and paste it in a web browser to get the same map?


texts: Raster Calculator and Land Change with multiple classifications

texts:
I am looking at land cover data that is separated into three classes: upland, wetland, and open water.

I have the data for varying years. I am interested in seeing how the land has changed over the years.

I used raster calculator and set one raster != to the other raster to show all the areas that changed through the year, but I am interested in a deeper analysis.

I would like to see if a wetland changed to an upland or open water and vice versa.

Is there a way to use the raster calculator with specific attributes in order to get my desired result? Or is there another way to do this?

I am using ArcGIS 10.2.2 for Desktop.


texts: Number of Max Instances Per Core in ArcGIS Server

texts:
ArcGIS Server 10.5 Windows

This question has been beat to death, but like others have said, documentation is hard to follow and described differently in places.  I understand all of the concepts and what is done, so I am hoping this is a simple question.

For a service, in some cases max instances can be set to n + 1 (n number of cores).

Does this mean for a 3 machine site with 4 cores on each machine the number is 12 + 1, or 4 + 1?

E.g. 4 cores each of three machines

4  x 3 + 1 = 13 max instances
or
4 + 1 = 5 max instances.


texts:
I am new to Python and coding generally, and I cannot figure out how to set up a debugger that will work with Arcpy. I am trying to use PyCharm, but it always gets caught up with the fact that it can't actually open ArcMap. It gives me this error when I run in the IDE:

  RuntimeError: Object: CreateObject cannot open map document.

I run the scripts from a Toolbox and they don't run correctly. All I know how to do is comment out some lines and narrow the problem down by re-running. There's got to be a better way. Right?


texts: QGIS is not working with Postgres Plus Advanced Server?

texts:
We have a problem to integrate GIS software with our GIS database which is populated on Postgres Plus Advanced Server 9.2 from EnterpriseDB. We are using QGIS to maintain and analyze our GIS repository, the software is working just fine when the database in the postgres community version, but when the database in ppas; QGIS return the following error while trying to fetch the list of GIS layers:

  Erroneous query:

 SELECT l.f_table_name,l.f_table_schema,l.f_geometry_column,upper(l.type),l.srid,c.relkind
 FROM geometry_columns l,pg_class c,pg_namespace n
 WHERE c.relname=l.f_table_name
   AND l.f_table_schema=n.nspname
   AND n.oid=c.relnamespace
   AND has_schema_privilege(n.nspname,'usage')
   AND has_table_privilege('"'||n.nspname||'"."'||c.relname||'"','select')
 ORDER BY n.nspname,c.relname,l.f_geometry_column

  returned 7 [ERROR:  schema "dbms_profiler" does not exist]

And the following is the answer from EnterpriseDB:

  We looked at query and it seems, issue is the way that query built. By
  default, PPAS stores package names in namespace table and so it is
  retrieving the names of packages as name of schemas hence query is
  failing. It works in PostgreSQL. Though, you are in postgres
  compatibility, we still retains some packages.

  So we would recommend you to contact QGIS vendor on the same as it is
  not tested with PPAS.


texts: How to organize/cluster nodes moving on the same direction

texts:
The Problem

I am working on a problem of clustering for mobile nodes.

Initially, I though of clustering the nodes based on their angle of movement (relative to the x-axis). However this situation would fail in case that we have two groups of nodes moving on the same direction but on two parallel tracks. Their angle of movement would be the same, as such my approach would place both groups on the same cluster, which would be wrong (see grey nodes in the fig).

Besides the angle of direction, also the (x,y,z) components of the nodes are known. Their distance from the origin can be computed as well.

The Question

What type of math/trigonometry/geometry calculations can be used to make a reasonable clustering scheme here?


texts: QGIS/GRASS: forecast slope variation in slope by a new road

texts:
For preliminary project purposes, I wonder if there is a way to forecast and calculate with QGIS and/or GRASS the variation of slope dued by new road, starting with terrain DEM and a project sketch of the road:

(translation: Piano di campagna pendenza 50/60% = original slope 50/60
              stradello larghezza 2,50m = road width).

Could be enough to made this task even in the maximum slope line.


texts: Single value of raster displays as transparent?

texts:
I have a GeoTiff raster that I'm trying to style, but a single value (1 of 8) refuses to take on any color. It always displays as transparent. Any idea of why this is / how to resolve it?

Working in QGIS 2.14/GDAL/GRASS.


texts: How to remove only a single point after using the Integrate tool

texts: I have a layer of 7000 point features that were obtained through asset extraction utilizing LiDAR. A small portion of these features are what we consider duplicates and need to be removed, however they do not have the exact same coordinates as they are with in a few feet of one another. I used the integrate tool to move each point to their relative middle, but now I need to remove only one of the two points and am having trouble doing so. I attempted to use the delete identical tool based on the feature's xy positions, but it did not remove either point as their coordinates were not exactly the same. How can I do this?


texts: How to represent seasonal variability in a single global NDVI map?

texts: I have twelve monthly MODIS datasets for a 1 year period (http://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD13A2_M_NDVI&year=2015) and wish to combine all the data into a single map highlighting average values but also representing seasonally high NDVI values in the Southern and Northern hemisphere summers. I have so far used the cell statistics tool in ArcGIS to create an annual average from the data and secondly an annual maximum (maximum values in each cell from all 12 datasets). I am not sure whether the average NDVI well represents high summer NDVI values, whilst maximum NDVI seems to over-represent them. Does anyone have any advice on how to best represent this data? I require just a single map as wish to produce a weighted overlay with a number of other raster datasets.


texts:
I'm hoping someone can help me here.

I have a line feature and what I am trying to achieve is to determine is in which direction the line ends.

I understand that this could be achieved with 'Feature Vertices to Points' to determine start and end nodes in ArcGIS Info, however I only have ArcEditor.

Are there any alternative ways to achieve this?

Thanks in advance for your help guys,

ps. I'm using ArcEditor 10.1 & I have QGIS.


texts: How to Select Layer by Location of Feature in PYQGIS?

texts:
I'm tried to select layer by location with processing tool in QGIS python.
its working Layer to layer but how to do with feature?
Code:
for area in AreaLayer.getFeatures():

                processing.run("native:selectbylocation", {'INPUT':OtherLayer,'PREDICATE':[0],'INTERSECT':area ,'METHOD':1})
                print (OtherLayer.selectedFeatureCount())


texts: Update Cursor not updating field correctly

texts:
I've got an issue with an update cursor not working correctly. The kicker is that when I dissolve by the field that isn't updating correctly, the output is actually correct. I have the same syntax on a different feature class, and it's working correctly. I've double checked field types - everything is good. I can't figure out why this is happening...

relevant code:

#populate DSA fields above threshold with current/correct assignment
arcpy.SetProgressorLabel("Populating DSA_Revised field with DSAs above the LOC: " +    str(Threshold))
with arcpy.da.UpdateCursor(OriginalSA,[DSA_Field,DSARevised_Field,LOC_Field],whereClause_UpdateCursor) as cursor:
    for row in cursor:
        row[1] = row[0] #Populate DSA_Revised field with DSAs that are above the user specified threshold
        cursor.updateRow(row)
        Assigned_List.append(row[0])#append entry to a list that will be used in later loop


texts: QGIS how to activate the checkboxes of the constractions in the attribute form by Python code and add expressions?

texts:
I tried using this type of line for the edit form , but its only for versions 2.8 but I don't know how is for new version 3.0 including "Enforce not null constraint".
formConfig = layer.editFormConfig()
formConfig.setNotNull( idx, True )


texts: Looping through multiple rasters, and adding them together

texts:
I have 5 folders, each with 80 .tif files. These 80 files represent one landsat image taken on a different date. Each folder represents a vegetative indice (NDVI, NDII etc.)

Each file in each folder is listed in the same order, based on date.  So the first file in all five folders correspond to each other.  I want to add all 5 .tif files together (for all 80 dates) using map algebra. The code I am trying so far is this:

import arcpy
from arcpy import *
arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True

arcpy.env.workspace=r'F:\Sheyenne\Normalized_Indices\NDVI'
NDVIrasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'F:\Sheyenne\Normalized_Indices\NDII'
NDIIrasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'F:\Sheyenne\Normalized_Indices\RGR'
RGRrasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'F:\Sheyenne\Normalized_Indices\SWIR32'
SWIR32rasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'F:\Sheyenne\Normalized_Indices\SATVI'
SATVIrasters=arcpy.ListRasters('*.tif*')
outpath='F:\Sheyenne\Normalized_Indices\Composited\\'
for ndviraster in NDVIrasters:
    ras1 = Raster(ndviraster)
    for ndiiraster in NDIIrasters:
        ras2=Raster(ndiiraster)
        for rgrraster in RGRrasters:
            ras3=Raster(rgrraster)
            for swir32raster in SWIR32rasters:
                ras4=Raster(swir32raster)
                for satviraster in SATVIrasters:
                    ras5=Raster(satviraster)
                    outRaster = ras1 + ras2 +ras3 +ras4 + ras5
                    outpath=outpath + ras1
                    outRaster.save(outpath)

but I think there is something wrong with how I am doing the nested loop, first of all, when I simply print ras5 at the end, it lists more than the 80 files found in the pathway r'F:\Sheyenne\Normalized_Indices\SATVI'

The actual error returned is:

Traceback (most recent call last):

  File "<ipython-input-61-567d86ebd99a>", line 1, in <module>
    runfile('F:/python codes/normalize_indices.py', wdir='F:/python codes')

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 685, in runfile
    execfile(filename, namespace)

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 71, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)

  File "F:/python codes/normalize_indices.py", line 98, in <module>
    outpath=outpath + ras1

  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 4143, in Plus
    in_raster_or_constant2)

  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)

  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 4140, in Wrapper
    ["Plus", in_raster_or_constant1, in_raster_or_constant2])

RuntimeError: ERROR 000732: Input Raster: Dataset F:\Sheyenne\Normalized_Indices\Composited\ does not exist or is not supported

but that pathway definitely does exists for the outpath.  Possibly an indenting issue in the last three lines?

EDIT:

I figured out the problem with the loop based on below answers and comments.

I am now using this:

outpath='F:\Sheyenne\Normalized_Indices\Composited\\'
for ndvi, ndii, rgr, swir32, satvi in zip(NDVIrasters, NDIIrasters,RGRrasters, SWIR32rasters, SATVIrasters):
    outraster= ndvi +ndii + rgr + swir32 +satvi
    outpath=outpath + ndvi
    outraster.save(outpath)
    print outpath

but the out path is using iteration to ad the raster names to the pathway. So the first outpath is the first raster in my list.  But the second one is the first and second combined, and the third one is the first through third combined etc.


texts: Is it possible to give 2d contours a Z-height from their attribute table in QGIS? (Or ArcGIS?)

texts:
I have a contour data file downloaded from UK Ordnance Survey, which gives 2d contours of an area in gml format (also supplied in .shp format) with the heights of each contour attached to the attributes of each line.

I want to give the contours a Z value of the height attribute attached to it, so that I have a 3d contour file I can then use in QGIS to create hillside shading etc, and potentially export to dxf/dwg.

Is there a way to do this without manually moving each line?

I was hoping for an actual tool I can use, as I have no script experience - I am willing to learn this if necessary, but would appreciate it if any answers that might refer to scripting please be as patient as possible, as I am clueless with where to even start to run a script!

I'm using QGIS at the minute but also have access to ArcGIS 9.2 (though none of the add-on program packages, just basic ArcMap.)

Thanks - Lisa


texts: Filling column with consecutive numbers in QGIS?

texts:
I make a new column in my attribute table and it has a default value (for example 0).
I want to (probably using the field calculator) have the result that row 1 has the value 1, row 2=2 and so on, just like a numbered index.

I'd be thankful for a code example for the field calculator (including the use of $rownum (I think I need this for switching to the line)).


texts:
I have no idea if this is even possible. We run a Utility software company and need to know where to put the leak sensors. The sensors are put at either every other house or every two houses depending on the size and material of the water pipe. We have all of the data (size of the pipes, material, etc.) in our GIS server and the meters and water lines are being plotted on a map.

Is it possible for us to be able to create an automated process for this?

Basically, one of our employees would go on our map and be able to click a button on our map and it would automatically plot these points based on the pipe material and the pipe size.

Can you plot every other address is ESRI?
Will I need a database of all the addresses in the United States?
Is there a basemap in ESRI with all of these addresses?
Or am I making this way harder than it actually is?


texts:
I'm working on a class project where we make a pie chart out of an attribute table.  The attributes are geologic layers (UNITNAME), and the area they occupy (area_pct).  In the "Create Graph Wizard", I have selected the correct table, and set the value field to "area_pct" and the label field to "UNITNAME".  However, I can't get the sort field to let me select anything.  I want the sort field to be "UNITNAME" as well.  I managed to get it to work yesterday by goofing around with the options, but I honestly have no idea how I did it, and now I can't replicate it.

Does anyone have any suggestions?  Is this arcGIS (10.2.1) just being picky as usual?


texts: Convert Sentinel GeoTIFF file to JPEG image

texts:
I am facing difficulty converting between Sentinel GeoTIFF file and JPEG file. I used below command i.e.

gdal_translate -b 4 -b 3 -b 2 -of JPEG -scale inputimage.tif outputfile.jpeg

The problem with this approach is, that if I convert 32*32 TIFF file into JPEG and if the same patch is from 1000*1000 TIFF file and I convert that file into JPEG, then in both cases the same patch looks different in both JPEG files.


texts: Is there a way to use a spatial query to invert X and Y coordinates in QGIS?

texts:
Original problem : I tend to often invert X & Y in WGS84 lat/long CSV import.

I'm wondering now if there is a way not to reimport the original file (as this is the basic solution for points) and invert geometries . My idea : a query on the geometry field may work to create a new layer based on the first one, but entering y,x coordinates instead of x,y in the geometry creation process, but I don't know yet how to achieve this.

Do you have an idea how I could use the powerful QGIS query engine to create new features from my original layer with inverted X & Y ? Answer can be specific for point/line/polygon as each of them may have a different solution type.

I'm using QGIS LTR (2.14) and latest (2.18)


texts: Using XSLT conversions in GeoNetwork to add XML nodes to metadata record during import

texts:
I want to use "Apply XSLT converison" option in GeoNetwork when importing new records to add new XML snippet to my metadata record.

I need to add aditional XML content (gmd:identifier) to the following node: MD_Metadata/identificationInfo/MD_DataIdentification/citation/CI_Citation

The XML that I want to add looks like this:

<gmd:identifier>
  <gmd:RS_Identifier>
    <gmd:code>
      <gco:CharacterString>0105</gco:CharacterString>
    </gmd:code>
    <gmd:codeSpace>
      <gco:CharacterString>hr:nipp:hr</gco:CharacterString>
    </gmd:codeSpace>
    <gmd:version>
      <gco:CharacterString>1.0</gco:CharacterString>
    </gmd:version>
  </gmd:RS_Identifier>
</gmd:identifier>

The MD_Metadata/identificationInfo/MD_DataIdentification/citation/CI_Citation node already contains one gmd:identifier tag, new one needs to go below it.

I need XSLT code that adds the code to specified node.


texts: Project data between different datums using custom parameters in QGIS!

texts:
I want to convert or reproject in the fly some shape files and raster between between different datums using custom transformation parameters (7 parameters) but don't know how to do that in QGIS.

I want to know the string of the final transformation, transformation matrice of rotation (my parameters are in seconds and scale is in ppm) that QGIS uses and if that string works only in one direction of the transformation or it automatically calculate the reverse transformation too?

In my scenario i want to go through this:

Epsg:3127 to Epsg:32634 (Transf. parameters are tx=0.0527 m ty=0.0509 m tz=-0.06636 m rx=-0.001456" ry=-0.008809" rz=0.014238" scale=0.000958 ppm)

Epsg:32634 to Epsg:28404 (Transf. parameters are tx=44.183 m ty=0.580 m tz=-38.489 m rx=2.3867" ry=2.7072" rz=-3.5196" scale=8.2703 ppm)

And how can i combine this 2 transformation paths in one so i can go directly from Epsg:3127 to Epsg:28404 and vice versa.

Here is the .gtf (32634 to 2642) contents created from the tool "create custom geographic transformation"

1DøÁÒ½  øu¼Î ]  GEOGTRAN["H_TRANS",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],GEOGCS["GCS_Pulkovo_1942",DATUM["D_Pulkovo_1942",SPHEROID["Krasovsky_1940",6378245.0,298.3]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],METHOD["Coordinate_Frame"],PARAMETER["X_Axis_Translation",44.183],PARAMETER["Y_Axis_Translation",0.58],PARAMETER["Z_Axis_Translation",38.489],PARAMETER["X_Axis_Rotation",2.3867],PARAMETER["Y_Axis_Rotation",2.7072],PARAMETER["Z_Axis_Rotation",-3.5196],PARAMETER["Scale_Difference",8.2703]]

and here is the .gtf (32634 to 3127 or 6870) contents created from the tool "create custom geographic transformation"

1DøÁÒ½  øu¼Î u  GEOGTRAN["UTM_KRGJSH_WGS_GK20FIN",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],METHOD["Coordinate_Frame"],PARAMETER["X_Axis_Translation",0.0527],PARAMETER["Y_Axis_Translation",0.0509],PARAMETER["Z_Axis_Translation",-0.06636],PARAMETER["X_Axis_Rotation",-0.001456],PARAMETER["Y_Axis_Rotation",-0.008809],PARAMETER["Z_Axis_Rotation",0.014238],PARAMETER["Scale_Difference",0.000958]]


texts: Troubles Updating Values of Selected Features with PyQGIS

texts:
I'm trying to  modify the values of some selected features. I have written this small script, but nothing happens when I run it. Is there something that I've missed? Here is my code.

lyr = qgis.utils.iface.activeLayer()
features = lyr.selectedFeatures()

lyr.startEditing()
lyr.updateFields()
for f in features:
    fid= f.id()
    lyr.changeAttributeValue(fid,12,"0")
    lyr.updateFeature(f)

print "elemento alterado"

lyr.commitChanges()


texts: Assigning point IDs to respective start and end attributes of polyline?

texts:
I have two shapefiles, Point and Polyline. I'm trying to populate my polyline layer with the point data ID (systemID) where the lines start to where they end; eg. create two columns (from and to). The lines already have flow direction assigned to them. I am using ArcHydro tools but Arcpy is also fine.

In ArcHydro I have tried using Attribute Tools > Generate From/To Nodes for lines but it creates values for To and From that are not from my point data. My other alternative was using Feature to Vertices and joining the result to my point data table to get the systemID. But then I would have to update the polyline layer as well. I am sure there is a easier way of doing it.

Note: not all the lines have points at the vertices, blank values are fine for them.


texts:
I have script that should add lyr file from server.
It is like a multiple symbology to certain layerin TOC, i.e. I have a feature class called tline and different symbology layers like tlineA, tlineB.
My goal is to have independent tline layer and upload tlineA and tlineB layers to TOC and change their data sources.
I made a script tool that already can read paths to tline and tlineA, tlineB and should change the data source of tlineA

Here is it

ab = arcpy.Describe("tline")
path = ab.path
ao = arcpy.Describe("tlineA")
peth=ao.path

for lyr in arcpy.mapping.ListLayers(mxd, "",df):
     if lyr.name == "tlineA":
         lyr.findAndReplaceWorkspacePath (peth,path)
         print lyr.dataSource

But anyway the datasource of tlineA remains the same as it was before.
I also tried like

lyr.findAndReplaceWorkspacePath (r"ao",r"ab")

But still doesn't work.


texts:
I have a digital elevation model produced from a drone survey via Drone Deploy. The file is a georeferenced TIFF.

The z values in the raster are relative to the survey only and not to actual ground level. So as an example, a z value of -39 on the raster is equivalent to 385m altitude in real life. 6.5 is 431m, -12 is 414m and so on. The difference is fairly consistent at 424 give or take.

I assume it would be a straight forward operation to apply a correction using the QGIS Raster Calculator.

What would the procedure be?


texts:
I am using pgr_dijkstra(or pgr_dijkstraCost) to find the shortest path and/or cost for drawing isochrone map. Since I need to calculate the shortest path from each of over 100,000 nodes to multiple destinations, I am planning to use pgr_dijkstraCost many-to-many function ↓

http://docs.pgrouting.org/2.2/en/src/dijkstra/doc/pgr_dijkstraCost.html

And, as shown in the example below, I have no problems when I input start_vids and end_vids as a form like ARRAY[2, 7] or ARRAY[3, 11].

SELECT * FROM pgr_dijkstraCost(
    'select id, source, target, cost, reverse_cost from edge_table',
    ARRAY[2, 7],
    ARRAY[3, 11]);

However, what I really want to do is to input arrays by selecting a column from a table. In other words, I want to calculate the shortest path from every single node in my table. For example, I've tried that

SELECT * FROM pgr_dijkstraCost(
    'select id, source, target, cost, reverse_cost from edge_table',
    select array(select source from edge_table),
    ARRAY[3, 11]);

But this returns an error: "ERROR: syntax error at or near "select" / SQL state: 42601 / Character: 171." What's wrong with my code? And, please tell me how to deal with the query result as an input parameter.


texts: Create a new Dataframe using ArcPy in 10.1

texts: Does anyone know if there is a way to create a new empty dataframe in ArcMap 10.1 using Arcpy?


texts:
In my php+postgis+geoserver+OL application I have to export a map in a pdf file.
This map has two layer: the base is a google map layer and a vector layer as overlay.
Is there a way to make a pdf like this?


texts:
On a map with territorial units, I wish to be able to compare a property of one territorial unit with those surrounding it.
How can I find the contiguous territorial units of a specific one ?


texts: Defining the geometric predicate as a dropdown-input of a model

texts: I built a model and use the tool "Extract by location". I would like to define the geometric predicate as a variable, so the user can choose between different predicates when he/she runs the model. How can I do that?


texts: Print progress using arcpy for blocking function

texts:
I am using ArcGIS 10.1 and through python I am calling arcpy's Clip function to clip a HUGE road network. This takes almost 2 days.

Is there a way I can see the progress?

Clip is a blocking call so I am not sure how to use SetProgressor.


texts: Maximum scale dependent visibility changes from 1:0 to -2147483648:1

texts:
QGIS 2.0.1 Dufour.
When i open the layer properties and the scale dependent visibility is unchecked and when i check it, the maximum (inclusive) shows 1:0. I click Apply, OK, reopen layer properties and it says -2147483648:1.
Everytime i change it to 1:0, apply and reopen the layer properties, the maximum (inclusive) gives me the value -2147483648:1.
Any ideas why?


texts: Determining optimal hexagon size for hexbins?

texts:
I've grown rather fond of hexbins. However, I'm always stumped when it comes to setting the hexagon size. I mostly go by feel right now. I just tweak the size until it "looks good".

Is there a more scientific way of determining the hexagon sizes?


texts:
Image corner coordinates are known:

Top left: (70.20315, 20.92749)
Top right: (70.20315, 31.45240)
Bottom left: (59.70571, 20.92749)
Bottom right: (59.70571, 31.45240)

(Basically map over Finland and a little bit of Sweden/Norway/Russia)

I am not sure what is current projection name; +1 in latitude and longitude is constant distance in image (which is not the case with Mercator projection).

I couldn't find a way to change image projection to Google maps overlay image, preferably using Python. Python isn't mandatory, but command line interface is.

It's possible to do that with pyproj and PIL, but I wasn't able to figure out how to do that. I also found GDAL, and more specifically gdalwarp, but couldn't find options for changing projection in a correct way.


texts: AERSURFACE - How to calculate the inverse distance weighted geometric mean of surface caracteristics

texts:
I need to calculate the inverse distance weighted geometric mean of the cell values of different sectors within the 1-km radius from a weather station, using QGIS. I'll try to make it clearly below with images,

Generally speaking, I have a weather station located on a random terrain (values of surface roughness 1, 2, or 3, for instance). For each section (0-90º, 90º-180, or other angles), within a radius of 1 km, I should calculate the distance-weighted inverse geometric mean. The average is calculated from the values associated with each cell of the land cover class within the individual sector considering the radius of 1 km from the meteorological tower. Below there is a rough example.

The values associated with each cell in the grid are weighted based on the inverse distance from the weather tower. Following the following formula:

I have no clue whether I can use the raster calculator to have this IDW Geometric mean for each sector, and how to do this, or if there is any other processing tool that could help me with?

Obs: there is a USEPA free software that calculates this, AERSURFACE, but it requires to the raster to be in Albers Conic Equal Area projection of USA, and when I try to use other Albers Conic Equal Area projection for areas outside the USA, then it doesn't run, and I think I cannot use the USA projection for non-USA areas, because I would have big deformations.


texts:
I found a way that should allow to use (Geo-)JSON flies with GeoServer: An unsupported GeoTools GeoJSON Datastore (thanks to this answer)
Here is a link to the jar: https://mvnrepository.com/artifact/org.geotools/gt-geojsondatastore/24.2
I dropped the jar into the WEB-INF => lib of GeoServer, but nothing happend to the GeoServer interface.
What else do I have to do to make the GeoJSON Datastore running to use JSON files with GeoServer?


texts: "Select features by expression" not matching all results for text field?

texts:
I have a vector layer containing all of the property parcels in the county. There is a field called "PROP_TYPE" that contains the type of property. I am trying to select all of the residential properties, which are coded as "RES" in the PROP_TYPE field.

I tried to do this by clearing current selection, making sure the parcels layer was the active layer, and then going to the "Select features by expression" dialog, and typed in the following queries:

"PROP_TYPE" IS "RES"

"PROP_TYPE" LIKE "%RES%"

The first returns only 533 features, but there are approximately 75,000 of them that should be selected. I can count them when sorted by prop_type in the attributes table, and see that they all contain exactly the same text in this field ("RES") ... so I can't figure out why some are matching and others aren't. Also, I don't understand why the first result would return 533 while the second would return 0. I would think the 2nd one would always have at least as many as the first, if not more, since the pattern includes more possible matches.

Can someone explain to me what is going on?