texts: Apply accuracy measures from partial area to the whole classified maps?

texts:
I applied supervised maximum likelihood classification to a yearly stack of Landsat images (resolution 30m). The land cover classes are related to coniferous forest, i.e., forest, clear-cut, fire, bark beetle. I had only few years of available reference aerial photos (resolution 0.5m).
Aerial photos do not fully cover the extent of Landsat images (AOI). To assess the accuracy of the classification for single year scene, I extracted a overlapping part between the extent of Landsat and aerial photos:

I calculated the measures of accuracy for the extracted part of fully classified Landsat by set of stratified sampling points per mapped class.

I followed the methods published here:

Olofsson, P., Foody, G.M., Herold, M., Stehman, S. V., Woodcock, C.E., Wulder, M.A., 2014. Good practices for estimating area and assessing accuracy of land change. Remote Sens. Environ. 148, 42–57. doi:10.1016/j.rse.2014.02.015 and
Olofsson, P., Foody, G.M., Stehman, S. V., Woodcock, C.E., 2013. Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation. Remote Sens. Environ. 129, 122–131. doi:10.1016/j.rse.2012.10.031

to assess the accuracy of the classification.

Great R code summarizing this can be found here:
https://github.com/openforis/accuracy-assessment/blob/master/error_matrix_analysis.R

The error matrix, by this method, is not based on the simple sample counts, but on the estimated proportions of the area. Including the proportion of the mapped classes in accuracy assessment, it is possible to quantify the uncertainty attributable to variability in the sampling. Simply, to adjust the mapped area to map classification error. The error-adjusted area estimates should be accompagnied by confidence interval.

My problem:

I classified 7 Landsat images. I calculated accuracy for extracted parts of them for two years, because of no availability of another reference data sets.

My question:

Can I apply the accuracy measures established from the partially covered area (intersection between extent of Landsat and reference aerial photography) to the whole classified scene Landsat?


texts:
Oftentimes when we make maps it is based on our subjective interpretation of what is aesthetically pleasing. I would like it if people posted examples of beautiful maps, displaying any phenomena in any manner.

Below I have posted one of my favorite maps. This is an example of a value-by-alpha map recently asked about as How to implement value-by-alpha map in GIS?, and the picture is taken from the GeoVista website.

Citation for the map's makers:

  Geovisual analytics to enhance spatial scan statistic interpretation:
  an analysis of U.S. cervical cancer mortality Jin Chen , Robert E Roth
  , Adam T Naito , Eugene J Lengerich and Alan M MacEachren
  International Journal of Health Geographics 2008, 7:57

It would be best for the cultivation of knowledge if people would elaborate on why the particular maps they cite are beautiful.

The reason I believe I think the cited value-by-alpha map is beautiful is that it creates a very simple, but obvious and striking visual hierarchy with which to interpret the standardized mortality ratio's. This is in particular useful combined with the very "noisy" standardized mortality ratio's, and the typically very noisy clusters of abnormally high rates produced by the SatScan clustering technique. One can even clearly see very small clusters around Chicago and Philadelphia.

There are other supplemental elements of the map that make it easy on the eyes. For example, the black background, the heavier white outline for around the states and the white outline for the states (that is blended the same as the attribute values). Maps with many polygons can particularly be distracting if one does not take care when plotting the polygon outlines.

Also the legend is particularly well created, and effectively demonstrates the concept (although it certainly isn't a typical legend, so took some original creative thought).


texts: Making feature layer selectable in ArcObjects Java without getting AutomationException 0x80004005?

texts:
I am currently making a FeatureLayer from a resultset.

The feature layer seems to work properly and will display on the map.

However when I choose "select all" on my feature layer I get the following automation exception:

AutomationException: 0x80004005 - Unspecified error
at com.esri.arcgis.carto.FeatureLayer.selectFeatures(Unknown Source)
at cc.nisc.mapping.controls.actions.selection.SelectAllAction.doActionPerformed(SelectAllAction.java:32)
at cc.nisc.gui.action.NiscAction.actionPerformed(NiscAction.java:225)
at javax.swing.AbstractButton.fireActionPerformed(AbstractButton.java:2018)
at javax.swing.AbstractButton$Handler.actionPerformed(AbstractButton.java:2341)
at javax.swing.DefaultButtonModel.fireActionPerformed(DefaultButtonModel.java:402)
at javax.swing.DefaultButtonModel.setPressed(DefaultButtonModel.java:259)
at javax.swing.AbstractButton.doClick(AbstractButton.java:376)
at javax.swing.plaf.basic.BasicMenuItemUI.doClick(BasicMenuItemUI.java:833)
at javax.swing.plaf.basic.BasicMenuItemUI$Handler.mouseReleased(BasicMenuItemUI.java:877)
at java.awt.AWTEventMulticaster.mouseReleased(AWTEventMulticaster.java:289)
at java.awt.Component.processMouseEvent(Component.java:6516)
at javax.swing.JComponent.processMouseEvent(JComponent.java:3320)
at java.awt.Component.processEvent(Component.java:6281)
at java.awt.Container.processEvent(Container.java:2229)
at java.awt.Component.dispatchEventImpl(Component.java:4872)
at java.awt.Container.dispatchEventImpl(Container.java:2287)
at java.awt.Component.dispatchEvent(Component.java:4698)
at java.awt.LightweightDispatcher.retargetMouseEvent(Container.java:4832)
at java.awt.LightweightDispatcher.processMouseEvent(Container.java:4492)
at java.awt.LightweightDispatcher.dispatchEvent(Container.java:4422)
at java.awt.Container.dispatchEventImpl(Container.java:2273)
at java.awt.Window.dispatchEventImpl(Window.java:2719)
at java.awt.Component.dispatchEvent(Component.java:4698)
at java.awt.EventQueue.dispatchEventImpl(EventQueue.java:735)
at java.awt.EventQueue.access$200(EventQueue.java:103)
at java.awt.EventQueue$3.run(EventQueue.java:694)
at java.awt.EventQueue$3.run(EventQueue.java:692)

my relevant code is listed below:

 try {
        String sqliteFile = String.valueOf(mwfFileSystem.getFesDirectory()) + "\\query_exports\\" + queryName + ".sqlite";

        workspace = WorkspaceFactory.SQLITE.openFromFile(
                sqliteFile);

        IFeatureLayer featureLayer = new FeatureLayer();
        featureLayer.setSelectable(true);
            featureLayer.setFeatureClassByRef(workspace.openFeatureClass(queryName));
        featureLayer.setName(queryName);
        ActiveMapContainer mapContainer = ActiveMapContainer.getInstance();

        mapContainer.getMap().addLayer(featureLayer);
    } catch (ArcEngineException | IOException e) {
        e.printStackTrace();
    }


texts: How to solve "Collection query aborted after accumulating over 5000 elements"?

texts:
I am trying to sample image and to get scatter plot of its' two bands- NDVI and db from senttinel 1.
For some reason, probably because of the size of the shapefile, I get the error:

  Error generating chart: Collection query aborted after accumulating
  over 5000 elements.

  Collection query aborted after accumulating over 5000 elements

That's happens after I try to sample my image using this:

// Generate a sample of points within the region
var sample = pairedImage.sampleRegions(geometry, null,3);

I have playes with the last number in the funnction , I started with 50 but then went down to 30,10,5 ,1 and 3, but non of them worked and in all I have gotten the same error. is it possible that it happens because of the size of the image and there is nothing I can do?

Here is a link to my code:

https://code.earthengine.google.com/f967f358a45cc1c40ccc6d0373066270

My end goal: to be able to sample this image and create the chart


texts: Calculating new attribute based on condition in QGIS 3 Field Calculator

texts:
I have a shape with different city areas that contains the average building heights ("HHO-mean") in meters. Now I want to categorize the average bulding heights in a new column in five different categories (1-5).
I tried the following:
CASE

    WHEN "HHO_mean" > 5.7 THEN 1
    WHEN "HHO_mean" > 6.4 THEN 2
    WHEN "HHO_mean" > 7.1 THEN 3
    WHEN "HHO_mean" > 8.0 THEN 4
    WHEN "HHO_mean" > 10 THEN 5

END

and
if( "HHO_mean" > 5.7, '1', if( "HHO_mean" < 6.4, '2',  if( "HHO_mean" < 7.1, '3', if( "HHO_mean" < 8.0, '4', if( "HHO_mean" < 10, '5', '0')))))

In both cases the new column only contains 1 in every field. Any suggestions?


texts: Thresholding Hansen Forest Cover/Global Forest watch treecover in earth Engine

texts:
I want to threshold the treecover2000 layer (i.e., tree cover in year 2000) at 25% using GEE. How would I do that? Below is code I wrote to load the gfc tree cover dataset

var gfcImage = ee.Image('UMD/hansen/global_forest_change_2015');
var treeCover = gfcImage.select(['treecover2000']);

Do I have to do something like below?

var treeCover25=gfcImage.select(['treecover2000']).select(0.25);


texts:
I have two datasets on fusion tables. One is a points layer (libraries) and the other is a polygon layer (land parcels). How do I design a query using Google Maps/fusion table to figure out this question: how many land parcels are within 1km (for example) of the library?

Also, I've seen some of the like this one where the search distance is pre-defined. Is it possible to design a query that doesn't have predefined search distances? I understand that this task will require the use of javascript (I'm still a beginner), so help in that area would also be greatly appreciated.

Thanks for your help!


texts: Automatically filling area surrounded by other polygons?

texts:
I want to fill the gap with a polygon automatically. It is surrounded by other polygons but the part inside the yellow circle is another shapefile. The gap is not made by adding a ring.

Or is my only choice to manually digitize?


texts:
I would like to scale polygon features by a specific factor I have calculated for each feature in a feature class.

I know of the Scale tool, however it only allows for manual entry of a scale factor and I have hundreds of features, so I would like to automate this as much as possible.

The end product will have all of my features resized so that they have the same area, but still have their own shapes and dimensions.

Does anyone know how I can accomplish this, preferably via python?  I can also utilize qGIS if a better function exists.


texts:
I have some problems for get and set properties in my features.

When I try get some property with

feature.get('property')

It not working. I have to look inside to the object and take the property like this:

feature["O"]["features"]["0"]["O"]["carrier"]

But it is not a good way, because when I changed CDN to archives in my app, the object is different. And, obviously, feature.get('') must work.

Also, when I try to set new properties to some feature.

var feature = map.getLayers().item(1).getSource().getFeatureById(id);
feature.setGeometry(new ol.geom.Point( ol.proj.transform([lon, lat], 'EPSG:4326', 'EPSG:3857') ));

It not working either.

Maybe the two mistakes are connected. And I suspect that it is for Cluster. It is the function when I create the Layer.

  var vectorSource = new ol.source.Vector({
      features: iconFeatures //add an array of features
  });

  var clusterSource = new ol.source.Cluster({
      source: vectorSource,
      distance: 40
  });

  var vectorLayer = new ol.layer.Vector({
       // source : vectorSource,
      source: clusterSource,
      style: clusterStyle
      });

    if (vector) {
      map.removeLayer(vector);
     }
     vector = vectorLayer;

      map.addLayer(vectorLayer);


texts: Can you run QGIS on MacOS High Sierra (10.13.3)? On the download page it says it's for "OS X Mavericks (10.9), Mountain Lion (10.8) and Lion (10.7)".


texts: Finding st_closestpoint between point to line using PostGIS?

texts:
I have point layer and line layer. I need to find the minimal distance between any point to any line.
I'm a beginner with postgres and I'm guessing I need to start with the st_closestpoint function. My question - Is there alternative of writing all the coordinates values?
the code example showing just a few of the coordinates, I Have to add many more to the function...

SELECT ST_AsText(ST_ClosestPoint(line,pt)) As closestpoint FROM (SELECT 'POINT(176684 660184)'::geometry As pt,
    'LINESTRING (178647 661077,  178722 661022, 178745 661006, 178760 661259, 178793 660970, 178795 661249, 178797 661196)'::geometry As line
) As foo;


texts: Converting Lambert Conformal Conic projected .las file to WGS84?

texts:
I am trying to import .las files in PostgreSQL database. The .las files are from https://viewer.nationalmap.gov/basic/?basemap=b1&category=ned,nedsrc&title=3DEP%20View. The issue is the file is in Lambert Conformal Conic projection and my remaining map data is in WGS84 projection.

I tried converting las file using las2las, but getting error:

  could not import coordinate system into OSRSpatialReference SetFromUserInput

I was following the tutorial https://gist.github.com/YKCzoli/3605e014b8ed09a571e5. But this tutorial does not talk about projection conversion.

Steps I followed:

las2txt -i SFdowntown.las -o SFdowntown.txt
copy elevation(x,y,z) FROM ‘~/Desktop/LiDARExtract/SFdowntown/SFdowntown.txt’ DELIMITERS ',';


texts: How would a Kriging model react to different z values at same point

texts:
Say I have a very large dataset of values collected at points over many years. For the sake of an example, we'll pretend it's the altitude at certain lat/longs.
If some of those points appeared twice in the dataset but with different z values (say due to a fairly extreme normal fault causing the altitude to change by a fairly large factor between one measurement and another several years later), then how would this affect an ordinary kriging model?


texts: Is there an API to ArcGIS Online's My Content listing?

texts:
This question pertains to ArcGIS Online, aka ArcGIS.com, aka ArcGIS for Organizations.

If I log into ArcGIS.com then choose My Content, I can choose > Layers > Features to see the feature layers available to my account.

Is it possible to get to this listing programatically, eg via a REST API to ArcGIS.com?


texts: QGIS Midvatten database error (python error while executing python code)

texts:
I am trying to create a Midvatten database and i keep getting error message that says "error has occured while executing python code (python error)".
Below is the error message:
2016-08-09T18:28:49 1   Traceback (most recent call last):
              File "C:/Users/agustus/.qgis2/python/plugins\midvatten\midvatten.py", line 654, in new_db
                newdbinstance = newdb(verno)
              File "C:\Users\agustus.qgis2\python\plugins\midvatten\tools\create_db.py", line 35, in init
                self.create_new_db(verno,user_select_CRS,EPSG_code)#CreateNewDB(verno)
              File "C:\Users\agustus.qgis2\python\plugins\midvatten\tools\create_db.py", line 78, in create_new_db
                self.rs = self.cur.execute(line.replace('CHANGETORELEVANTEPSGID',str(EPSGID)).replace('CHANGETOPLUGINVERSION',str(verno)).replace('CHANGETOQGISVERSION',str(qgisverno)).replace('CHANGETOSPLITEVERSION',str(versionstext[0][0]))) # use tags to find and replace SRID and versioning info
            UnicodeEncodeError: 'ascii' codec can't encode character u'\xf8' in position 8: ordinal not in range(128)

2016-08-09T18:45:19 1   Traceback (most recent call last):
              File "C:/Users/agustus/.qgis2/python/plugins\midvatten\midvatten.py", line 654, in new_db
                newdbinstance = newdb(verno)
              File "C:\Users\agustus.qgis2\python\plugins\midvatten\tools\create_db.py", line 35, in init
                self.create_new_db(verno,user_select_CRS,EPSG_code)#CreateNewDB(verno)
              File "C:\Users\agustus.qgis2\python\plugins\midvatten\tools\create_db.py", line 78, in create_new_db
                self.rs = self.cur.execute(line.replace('CHANGETORELEVANTEPSGID',str(EPSGID)).replace('CHANGETOPLUGINVERSION',str(verno)).replace('CHANGETOQGISVERSION',str(qgisverno)).replace('CHANGETOSPLITEVERSION',str(versionstext[0][0]))) # use tags to find and replace SRID and versioning info
            UnicodeEncodeError: 'ascii' codec can't encode character u'\xf8' in position 8: ordinal not in range(128)

i do have no idea how to fix this.


texts: Working with AIG Raster format in ArcGIS for Desktop?

texts:
I downloaded a specific layer which has the format "AIG" as you can see in the screenshot.

It has no projection so I applied "Define Projection" on it. I searched for the projection and found the following:

  4.1 - INDIRECT SPATIAL REFERENCE SYSTEM:

  Type of indirect spatial reference system:

  By country.

  Reference date:
   2001.

  4.2 - DIRECT SPATIAL REFERENCE SYSTEM:

  Datum:  Lambert Azimuthal.

  Ellipsoid: Sphere.

  Map projection: Projection    LAMBERT_AZIMUTHAL

  Units                    METERS
             Spheroid         SPHERE

  Major Axis        6378388.00000
              Minor Axis       0.00000

  Parameters:

  radius of the sphere of reference                      6378388.00000

  longitude of centre of projection                        9  0  0.000

  latitude of centre of projection                        48  0  0.000

  false easting (meters)                                       0.00000

  false northing (meters)                                      0.00000

The problem I have is that if search for this projection I get these two

After applying the first one using "Define Projection", nothing changes in the file. By this I mean the spatial reference is still "undefined" and the coordinate system is unknown. I also tried other projections but nothing works. If I convert it to a polygon first, define projection works but the reference on the map is still completely wrong. I think these problems occur because of the "uncommon" raster format?

I also think that the projection from the screenshot is not equal to the one I posted up here?

My goal is to project it in the end to "British National Grid" as all my layers I work with have this projection.

I also tried it on another computer with ArcGIs 10.2 and the same problem occured.

See these FME links for further reference to the AIG format:

Esri ArcGIS Binary Grid (AIG) Reader
Esri ArcGIS Binary Grid Quick Facts

I am using ArcGIS 10.2.2 for Desktop.


texts:
Can anybody help me to solve a problem with running gdal_calc for many rasters?

I have several multispectral images and I need a batch to perform NDVI calculation.

I started with gdal code for one raster and it works good:

> C:\Python27\python.exe "C:\Program Files\QGIS Wien\bin\gdal_calc.py" -A D:\my_dir\ras0_2.TIF --A_band=3 -B D:\my_dir\ras0_2.TIF --B_band=4 --outfile=D:\my_dir\ndvi.tif --calc=((B-A)/(B+A))

But when I tried to run code for all TiFFs in a folder the program gave an error
"nonetype object has no attribute GetRasterBand"

for %i in (*.tif) do C:\Python27\python.exe "C:\Program Files\QGIS Wien\bin\gdal_calc.py" -A $i --A_band=3 -B $i --B_band=4 --calc="((B-A)/(B+A))" %i ndvi/ndvi_%i


texts: Is it possible to create temp or infile variable instead of creating a new feature?

texts:
I am creating a small script which requires to run a Buffer_analyst tool. The help of this function says that I have to set a path where the function will save the output. Is it possible to store the output in a file/script/fuction variable or a temp file instead of saving it permanently to my drive?

The help says that I have to do this:

arcpy.Buffer_analysis(input_file, output_file, buffer)

but I would like to get something like this:

output_file = arcpy.Buffer_analysis(input_file, buffer)

or this:

arcpy.Buffer_analysis(input_file, temp_file, buffer)


texts: Repairing all topology errors at once in ArcMap?

texts:
In 2012 an ArcMap user wrote Automatically fixing multiple errors using Fix topology error tool?:

  I am a ArcGIS.10. user. I'm want to solve multiple errors using the
  fix topology error. So I select all my features I click right but I
  don't have the correcting options (Snap, trim, extent), I only have
  Delete. Does anyone know why?

Now I have about the same problem with overlapping line features which I would like to correct with subtract without the need to select and subtract each single error.

Is there in 2018 the option to correct all errors in one go?

For overlapping line features there is only the option to subtract, but still when I select multiple errors the option to use subtract is greyed out.

Currently I have more than 10000 errors in a big data set of an indoor environment.


texts: How to update circle-color or marker symbol on update of a source property in Mapbox GL?

texts:
I am trying to get Mapbox GL work in my React code as below:

I have a property called strength in each of my features in a FeatureCollection of a mapbox-gl map. I want to display a circle (ideally a marker) based on the value of strength property (0-3: green, 4-7: yellow, 8+: red).
To test this, I am initially having three features with strength as 3(green), 6(yellow), and 9(red) respectively. I have a button on click of which the strength of all the three features updates to 8. At this point, I want all the circle colors to become red. However, there is no change in the colors.

The click handler of the button is as below:

handleClick = () => {
  const newFeatures = this.state.features.map(feature => ({...feature}));
  newFeatures.forEach(feature => feature.properties.strength = 8);
  this.setState({feature: newFeatures}, () => this.state.features.forEach(feature => console.log(feature.properties.strength))
  );
}

The features are defined in the component state and componentDidUpdate looks like this:

componentDidUpdate() {
  console.log('Component Updated');
  this.state.features.forEach(feature => console.log(feature.properties.strength));
  this.map.getSource('institutions')._data.features.forEach(feature => console.log(feature.properties.strength));
  this.map.setPaintProperty("institutions-layer", 'circle-color', generatePropertyFunction());

generatePropertyFunction looks like this:

function generatePropertyFunction() {
  const propertyFunction = {
    property: "strength",
    type: "categorical",
    default: "rgba(0,0,0,0)",
    stops: [
      [0, '#7ff000'],
      [1, '#7ff000'],
      [2, '#7ff000'],
      [3, '#7ff000'],
      [4, '#ffff00'],
      [5, '#ffff00'],
      [6, '#ffff00'],
      [7, '#ffff00'],
      [8, '#ff0000'],
      [9, '#ff0000'],
      [10, '#ff0000']
    ]
  }
  return propertyFunction
}

This works the first time the map is rendered as the 3 locations are marked in green, yellow and red as intended. But after click of button, when the property strength is updated to 8 (this is verified in both console.log statements) the colors of the marker don't change. However, on click of button, if I explicitly set the circle-color to any one color, it works as below:

this.map.setPaintProperty("institutions-layer", 'circle-color', '#ff0000');

Am I missing something here or is this a bug?


texts:
# Declare instance attributes
self.actions = []
self.menu = self.tr(u'&City_searcher_plugin')
# TODO: We are going to let the user set this up in a future iteration
self.toolbar = self.iface.addToolBar(u'city_searcher_class')
self.toolbar.setObjectName(u'city_searcher_class')
QObject.connect(self.dlg.Statecombo, SIGNAL("currentIndexChanged(int)"), self.findTheCityNames)
QObject.connect(self.dlg.Citycombo, SIGNAL("currentIndexChanged(int)"), self.calcDistance)
QObject.connect(self.button_box, QtCore.SIGNAL(_fromUtf8("accepted()")), self.buffer)
#QObject.connect(self.dlg.browsepushButton, SIGNAL("clicked()"), self.saveTheFile)

In the above snippet, I used fromUtf8. But it is saying an error fromUtf8 is not defined. Where should I import fromUtf8 like from qtgui or qtcore?


texts: I have be finding a clear differences between trip and route terminologies. I am also developing a transit system using goole transit. This thread answers has a very clear explanation. But one confusion remains, can trip exist without routes? I think logically it can. But what about in google transit context?


texts: Find same names (duplicates) in one field using another feature class field

texts:
I have 2 tables (UML, PalmOil) and I would like to check record by record within PalmOil, if there are duplicates intersecting NameUMLwith NamePalmOil
UML
ID | NameUML |
1    NAME1
2    NAME2
3    NAME5

PalmOil
ID | NamePalmOil |
1    Name1
2    Name2
3    Name4

# Create a field to store all the PalmOil names in uppercase
arcpy.AddField_management(PalmOil, 'upper_name', 'TEXT')
# Calculate field
arcpy.CalculateField_management(PalmOil, 'upper_name', '!NamePalmOil!.upper()')

NameUML = 'NameUML'
UpperName = 'upper_name'

# Create an empty list to store all the Uml Names
UmlNames = []

#Create an empty list to store all the PalmOil Names
PalmOilNames = []

# fill in UmlNames list using a search cursor
fields = [NameUML]
with arcpy.da.SearchCursor(UML, fields) as cursor:
    for row in cursor:
        UmlNames.append(row[0])
del cursor

# fill in PalmOil name list using a search cursor
UpperName = 'upper_name'
fields = [UpperName]
with arcpy.da.SearchCursor(PalmOil, fields) as cursor:
    for row in cursor:
        PalmOilNames.append(row[0])
del cursor

# find identical values within the two lists
result = []
IdenticalNames = set(UmlNames).intersection(set(PalmOilNames))
result.append(IdenticalNames)
print (result)

So far so good:here if I print the list I can see all the matching names as I wish
Then I add a field to store the risk of duplicates, and I run an update cursor to update it.

# Once i have the results for same names i can upload the PalmOil dataset

# Add field 'DuplicateName'
arcpy.AddField_management(PalmOil, 'duplicate_name', 'text')

# Run an update cursor to update the duplicate name fields
# (this step is not working)
for i in result:
    print (type(i))
    fields = ['upper_name', 'duplicate_name']
    with arcpy.da.UpdateCursor(PalmOIl, fields) as cursor:
        for row in cursor:
            print (type(row[0]))
            if row[0] == i:
                row[1] = 'yes'
                cursor.updateRow(row)
            if row[0] != i:
                row[1] = 'no'
                cursor.updateRow(row)
    del cursor

The last cursor is working, but not as I want.
it just update the table with 'no' also if some of the name are matching according to the intersection I did before.
Where am I wrong?
This is the result i would like to have in PalmOil table:
PalmOil
ID | NamePalmOil | upper_name | duplicate_name |
1    Name1         NAME1         yes
2    Name2         NAME2         yes
3    Name4         NAME4         no


texts: Enabling toggle editing when loading layers form postgis in QGIS

texts:
I'm trying to load layers from PostGIS using the code bellow:
uri = QgsDataSourceUri()
uri.setConnection(server, PORT, base, user, mdp)
root = QgsProject.instance().layerTreeRoot()
group = root.addGroup(groupname)
for element in table_list:
if element == "bpe":
    asql = """SELECT id,
                 geom,
                 nom,   # primary Key
                 code,
                 adresse,
                 cp,
                 ville,
                 type_fonc,
                 etat,
                 modele,
                 reference,
                 Proprietai,
                 gestionnai,
                 date_de_cr,
                 emprise
           FROM  public.bpe
           INNER JOIN public.bpe_projet AS bp ON bpe.nom = bp.nom_bpe
           WHERE bp.nom_projet = '%s'""" % project
    uri.setDataSource('', '(' + asql + ')', 'geom','','nom')
    layer = QgsVectorLayer(uri.uri(), element, "postgres")
    QgsProject.instance().addMapLayer(layer, False)

When the layers loaded into QGIS, the toggle editing is disabled.
While I'm using the code below, the toggle editing is enabled.
elif element == "support":
    asql = "nom_projet= '%s' " % project
    uri.setDataSource('public', element, GEOM, asql)
    layer = QgsVectorLayer(uri.uri(), element, "postgres")
    QgsProject.instance().addMapLayer(layer, False)

What modification I need to do for the first code for enabling the toggle editing?


texts:
I'm trying to rasterize the "Decreasing land productivity" shapefile from https://wad.jrc.ec.europa.eu/geoportal, using:

gdal_rasterize -3d -co COMPRESS=DEFLATE -ts 43200 21600 ./lpd_int2.shp ./ldp_int2.tif

(also tried without the COMPRESS=DEFLATE)
It has been running for 6 days and is only 40% completed.

Though the output is 43200x21600, this still seems like an excessively slow process. Is there some obvious thing I'm doing which is making it take too long?

$ ogrinfo -ro -al -listmdd lpd_int2.shp
INFO: Open of `lpd_int2.shp'
      using driver `ESRI Shapefile' successful.

Layer name: lpd_int2
Metadata domains:
  (default)
Metadata:
  DBF_DATE_LAST_UPDATE=1919-09-29
Geometry: Polygon
Feature Count: 226
Extent: (-179.999989, -55.600000) - (179.991667, 72.400000)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433],
    AXIS["Longitude",EAST],
    AXIS["Latitude",NORTH]]
Data axis to CRS axis mapping: 1,2
id_country: Real (33.15)
name_iso31: String (254.0)
iso3: String (254.0)
sqkm: Real (33.15)
lpd_sqkm: Real (33.15)
lpd_p: Real (33.15)
shape_leng: Real (33.15)
shape_area: Real (33.15)
OGRFeature(lpd_int2):0
  id_country (Real) = 3.000000000000000
  name_iso31 (String) = Afghanistan
  iso3 (String) = AFG
  sqkm (Real) = 643490.775693999952637
  lpd_sqkm (Real) = 52103.569999999999709
  lpd_p (Real) = 8.100000000000000
  shape_leng (Real) = 90294575.479599997401237
  shape_area (Real) = 80173734283.800003051757812
  MULTIPOLYGON  ...etc, etc...


texts: Where is the Production Editing Toolbar?

texts: I'm having trouble finding the production editing toolbar in arcgis 10.1. Is it part of an add on extension?


texts: Live update (the database & Map) possible in ArcGIS Server 10.1?

texts:
I am using SQL Server server database contains lat/lon values and using this information I have created the shapefiles and published on server...

When I get the update say 100 records then I need to recreate the shapefile and merge it to (previous one) like this I am doing daily for last few days...

It is difficult to do daily for all records...

So is there any workaround?? I am using ArcGIS Desktop 10.1 and ArcGIS Server 10.1+ MS sql database 2008 r2

Any help will be great!


texts: Changing attribute and/or text values in Esri metadata using Python

texts:
I'm new both to this forum and Python, and I'm having trouble tracking down a solution to my problem. I'm hoping someone can offer some process-improvement advice and/or code to help get my metadata issue sorted. This is the first script I've attempted so I imagine there's a better way to do what I'm looking for.

My Goal:

Create a script to:

•   import metadata and have it reflect the correct information about the new file
•   remove linkage, lineage, and copyHistory elements or text (and any other reference to location on our servers) for public sharing of file

Current process:

•   Python script to apply metadata based on an existing XML file using arcpy.ImportMetadata_conversion, sync enabled.
•   Remove elements from XML file using minidom removeChild.

# Step 4: Import metadata

# polygon file
arcpy.ImportMetadata_conversion ("metadatafile.xml", "FROM_ARCGIS", polyfilefull, "ENABLED")

#Step 5: Remove identifying information from metadata
# based on http://stackoverflow.com/questions/5699745/how-to-comment-out-an-xml-element-using-minidom-dom-implementation
# added "utf-8" to avoid a whitespace error, appears to work

from xml.dom import minidom

polydoc = minidom.parse(polyfilefull + ".xml")
    for linkage in polydoc.getElementsByTagName('linkage'):
    linkage.parentNode.removeChild(linkage)
for lineage in polydoc.getElementsByTagName('lineage'):
    lineage.parentNode.removeChild(lineage)
for copyHistory in polydoc.getElementsByTagName('copyHistory'):
    copyHistory.parentNode.removeChild(copyHistory)
f = open(polyfilefull + ".xml", "w")
f.write(polydoc.toxml("utf-8"))
f.close()

Issues:

•   “Linkage” reappears when metadata is viewed in Arc.
•   Want to update publication date based on current date

Possible solutions:

•   Remove the text of the element rather than the element itself, and remove ‘sync’ from the element or
•   change the value of sync to FALSE.


texts: Batch Clip in ModelBuilder of ArcGIS for Desktop?

texts:
I have multiple rasters and feature classes that I want to clip to a study area.

I am pretty new to ModelBuilder.

I know you can right click the clip tool to batch clip outside of ModelBuilder but how do I bring that into ModelBuilder at ArcGIS 10.1 for Desktop?


texts:
I am trying to find a windows installer for GeoServer.
Could not find it :(

Does anyone know where to find it or how to set up GeoServer on Windows using Platform independent binaries or is there another way to set it up?

I tried going to this page but could not locate windows installer.
http://geoserver.org/release/maintain/

I have used GeoServer earlier but at that time I was able to find a windows installer. Not sure if GeoServer has stopped providing windows installer.


texts: Replace NA's with 0 for large raster data using R?

texts:
I want to replace all NA values in a large raster with 0. I found this answer:
Replace NA's with 0 for raster data using R?

#getting a raster
library(raster)
f <- system.file("external/test.grd", package="raster")
f
r <- raster(f) #r is the object of class 'raster'.

# replacing NA's by zero
r[is.na(r[])] <- 0

but it's not working for my case because the raster is too large thus the computer runs out of memory. Is there an other solution which is more memory efficient?


texts: Hiding specific points in shapefile layer

texts:
I am using QGIS 2.18.2 on a Macbook Pro running OS X 10.10.5.

I created a layer of only 8 points from a wildlife survey. A few points have become less relevant, so I want to hide them from the display. But I would rather not delete them from the shapefile.

Is there a way to do this? Or do I have to create a new, separate layer of the relevant points?


texts: Using C# with ArcGISRuntime, how do I use a StretchRenderer to make a grayscale color ramp?

texts:
Using C# with ArcGISRuntime I need to display raster data on a map.  Right now I have the following code that successfully displays an elevation raster (DEM) on the map.  However, I have other rasters that are not displayed correctly using this color ramp.  When I bring them into Arc Map, they are displayed the way I want them, as grayscale images.  How can I define a StretchRenderer that would be a grayscale color ramp, that would make these rasters appear the same as they do when bringing them into Arc Map?

private void AddRasterLayer(string rasterPath, string id)
{
    var newRaster = new Raster(rasterPath);
    var newRasterLayer = new RasterLayer(newRaster);
    newRasterLayer.Id = id;
    // Code to get the min and max values for the MinMaxStretchParameters object
    double[,] rasterArray = Conversions.RasterToArray(rasterPath);
    double max = SoilMoisturePreprocessing.Max2d(rasterArray);
    double[] minVals = { 0 };// 0
    double[] maxVals = { (int)max + 1 };
    MinMaxStretchParameters minMaxParams = new MinMaxStretchParameters(minVals, maxVals);
    ColorRamp elevationColorRamp = ColorRamp.Create(PresetColorRampType.Elevation);
    var rendrr = new StretchRenderer(minMaxParams, null, true, elevationColorRamp);

    newRasterLayer.Renderer = rendrr;
    newRasterLayer.Opacity = 0.6;
    var mapViewModel = ServiceLocator.Current.GetInstance<MapViewModel>();
    mapViewModel.Map.OperationalLayers.Add(newRasterLayer);
}


texts: Road snapping in R - snapPointsToLine / Projection Problem

texts:
1. Introduction (you may skip this if you like)

In order to implement a road snapping feature in R for recorded positioning data, I first searched the internet for possible solutions to the problem:

While it seems this can be solved with fuzzyMM (using fuzzy logic), the project using that approach is unfortunately not active any more and, thus, the documentation does not include the kind of information that would be required to get started with the package itself.

Another solution could be using online geospatial APIs, e.g. Google's. However, the current database already includes way more datasets than what's included in the free package.

At last, I found maptool's function snapPointsToLines (short introduction about how the function works can be found here) and attempted to snap a couple of coordinates to lines provided by a shapefile.

Note: Because the GPS-data was recorded on cars mostly driving on major roads in Germany, I downloaded the respective shapefile (of Bavaria, i.e. Bayern) here (.shp.zip).

2. The specific problem

Road snapping with snapPointsToLines() generates quite confusing values with this code. The following picture shows a list of recorded longitude / latitude values (see pointsCoord).

roadsBavaria = readOGR("~/roads.shp", "roads")
pointsCoord  = getValuesFromPicture()

pointsSpatial          = SpatialPoints(pointsCoord)
pointsSpatialCorrected = snapPointsToLines(pointsSpatial, roadsBavaria)

Although this reduced test case sets the value pointsSpatialCorrected (see result at the end), it fails (code after that is not executed) with

Warning messages:
1: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
  Spatial object 1 is not projected; GEOS expects planar coordinates
2: In RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance") :
  Spatial object 2 is not projected; GEOS expects planar coordinates

Note that this error occurs even when I apply the shapefile's projection on the spatialPoints manually via

pointsSpatial = SpatialPoints(pointsCoord, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

What am I doing wrong here?

3. Additional information

The resulting value of spatialPointsCorrected as displayed by View()


texts:
Leaflet-routing: http://www.liedman.net/leaflet-routing-machine/

If you initialize the L.Routing.control like this:

var router = new L.Routing.control({
    waypoints: [L.latLng(57.74, 11.94), L.latLng(57.6792, 11.949)],
    language: 'pt' // here's the magic
}).addTo(map);

It appears translated on the page.

but when I access the instructions in the 'router' object: router._routes[0].instructions[0] it returns something like this:

router._routes[0].instructions[0] = {
    direction: "E",
    distance: 13041.6,
    exit: undefined,
    index: 0,
    mode: "driving",
    modifier: undefined,
    road: "",
    time: 860.2,
    type: "Head"
}

How can I get the translation formatted like in the interface?


texts: How to sum values for any x,y coordinate as key in a dictionary by PyQGIS?


texts: Using .len() in IF THEN statement in ArcMap field calculator

texts:
I have a string field with values of varying lengths (1-4).  I need each value to have a length of 5 characters by adding zeros to the value until the length equals 5.  For example:

123

needs to be:

00123

I am writing a Python function in ArcMap's Field Calculator but it is incorrect.  This is my code:

def zeros(photoLen):
  if photoLen == 4:
    return "0"
  elif photoLen == 3:
    return "00"
  elif photoLen == 2:
    return "000"
  elif photoLen == 1:
    return "0000"
  else:
    return 0

And then I call the function like this:

zeros(!photoStr!)

What do I need to change so this function will work?


texts: What is the schema of HEC-RAS output in HDF5 format?

texts:
I realize that HDF5 is a self describing format, and I can use a tool or programmatically find out what types of data groups exist in an HDF5 file.

What I am really looking for are the group and data set names created by HEC-RAS in the pxx and hxx file types. This information does not appear to be documented on the HEC-RAS web site or in the numerous manuals the group publishes.


texts: ArcGIS 10.2 Extract Polyline from Start/End Coordinates

texts:
I have two files. A polyline shapefile of highways and a spreadsheet that contains beginning and ending coordinates with descriptions of segments of the highways. For example:

BegX, BegY, EndX, EndY, DescriptionField1, DescriptionField2
-80.0, 30.0, -80.1, 30.1, Hwy1Partial, 1
-80.1, 30.1, -80.2, 30.2, Hyw1Partial, 2

I need to create a new polyline shapefile that breaks up the highway segments from the original polyline file into a new one that has the beginning and ending coordinates from the spreadsheet.

I've attempted to join the data but can't figure out how to get it to extract the polyline that starts at the begin coordinate and ends at the end coordinate. I've used the XY to Line tool, but that just gives me a line connecting the begin coordinate to the end coordinate. I need to match it up with the polyline from the highway shapefile.

The shapefile I am using does not have M values. The spreadsheet I have contain latitude/longitude coordinates that look to have originally came from the shapefile (they fall directly on top of it when plotted). I need a polyline file that matches the spreadsheet exactly. Meaning the spreadsheet would become the new shapefiles attribute table, and each record would represent a polyline that was extracted between the beginning coordinate and end coordinate.

The black line represents the shapefile that I have. The two coordinates are examples of beginning and ending lat/lon coordinates that are in the spreadsheet. What I want is to extract the blue line from the shapefile into a new polyline file with the attributes of the original spreadsheet.

I have played around with making XY point layers out of my beginning and ending coordinates then splitting the file based on those new points. The problem is that my shapefile has more lines than I actually need. The resulting file doesn't split it completely the way I need, and I cannot find a way to relate my attributes from the original spreadsheet back to it. What I really need is a custom query that takes in a begin lat/lon, then an end lat/lon and queries the shapefile to get the resulting line. After it finds the new line it would take it to a new shapefile; copying the attributes over from the spreadsheet to be associated with the newly created line.


texts: View and edit roles in Portal for ArcGIS and ArcGIS Web AppBuilder?

texts: I am using Portal for ArcGIS 10.4.1 and web app builder developer edition to build application. My application is ready but now I have to create two groups: One group have a view capability only and other can have edit capabilities in the same application.


texts: Seeking world map (epsg 3857) tile boundaries in shapefile?

texts:
I am working on manipulating data and I need a shapefile of world map tile in each zoom level.

It could be better if it is in polygon but line is also fine.
I tried to find the file using as many keywords as I can. Yet, I was not able to find anything.

It has to be as same as mapbox tile boundaries.
For those who have this data, please do share for the GIS users.


texts: Is there a way to quickly flip through PostGIS records with several geometry fields?

texts:
I often find myself in the situation of needing to look at several spatial attributes for a particular item at once, each associated with a single ID field or record, and wanting to look at several such records in sequence. For example, I might have some data where I have a set of say, a hundred users, each with a home location, work location, a GPS trace, and so on. Let's say these are all in the same row in a PostGIS table, but since they are separate geometry fields, need to be loaded in QGIS as several distinct layers. If I want to look at all the data for user 1, what I've been doing is setting an identical "provider feature filter" for each layer so I can filter the data to just that user. That's fine for one user, but if I then want to look at user 2, user 3 and so on I have to (for each layer) open up the layer properties window, open the query builder, and increment the id. This can get very tedious if each user has many geometry features!

The way I've solved this before is to create a new field in my PostGIS table(s) to indicate an 'active' user, setting those records to true and all others to false and then targeting that field with the layer filters. But I still have to flip to a new window and run a query after editing some text for each new record.

I'm wondering if there might be a simpler/quicker solution to this, without having to leave the QGIS window? Is there some way to link layers so they share a filter?


texts: ERDAS file not recognized as ERDAS Hierarchal File Architecture (EFHA)?

texts:
I am trying to do an accuracy assessment in Erdas using "user defined ground control points" (I collected the GPS points myselfin the field using a handheld GPS).

However, when I load the thematic image I want to validate, I get an error message saying the file is not an EFHA file (ERDAS Hierarchal File Architecture) which results in the file not being loaded into the CellArray.

Any idea on how I can I solve this problem?


texts: Count how many new points are captured when a polygon buffers by x miles

texts:
I have a polygon and a layer of many points. How many more points would be captured if I expanded the polygon by x meters? I want to end up with values for x = 1 mile, 5 miles, 10 miles, and maybe a few more. My instinct is to plug it into GIS, calculate a buffer, cut out the original polygon, and count point features, but it may not be feasible when I get more points: I will eventually be using a layer of tens of millions of points.
Is there a faster way to do this, maybe outside of a GIS?


texts:
Just wondering if there is a universal map format that all GIS software can read.
ie to quickly share layers with all our partner agencies which have different GIS Software, rather than have to convert it in multiple different formats (from our end or their end).
We use ESRI, but have partners using MapInfo, Blue8, CadCorp...

I wondering If I was missing a trick.


texts: Sinusoidal to WGS84 using Project Raster in ArcGIS

texts: I am planning to use the project raster tool in ArcGIS 10.7 to convert a custom sinusoidal projection (MOD10A1) to GCS_WGS_1984. I know I can use HEG or gdal (not familiar with gdal for now), but can I instead use the Project Raster tool for this type of coordinate conversion? If yes, then what kind of resampling technique is more accurate? I have attached a screenshot for further reference. The purpose of this transformation to align to the pixels of two different rasters, which currently have different CRS. My project area is in the Pacific North West (US). Basically, I want to confirm whether project raster is the correct method for this kind of coordinate conversion or not.


texts:
Possible Duplicate:
  Problems with OpenLayers and WFS - proxy.cgi issues

I have a problem with a simple OpenLayers web page. I have two layers, first one is a WMS layer and work perfectly on QGIS and OpenLayers but the second WFS layer (tasmania_road) works on QGIS but it doesn't work on OpenLayers. I simply show two layer with  WFS tasmania_road on top of the WMS layer. My code shown below.

<html>
<head>
    <script type="text/javascript">
        var map;
        function init() {
            OpenLayers.ProxyHost = "proxy.cgi?url=";
            map = new OpenLayers.Map('rcp1_map');
            var wms = new OpenLayers.Layer.WMS(
                    "OpenLayers WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0",
                    { layers: 'basic' }
                );

            var layer = new OpenLayers.Layer.Vector("WFS", {
                strategies: [new OpenLayers.Strategy.BBOX()],
                protocol: new OpenLayers.Protocol.WFS({
                    url: "http://demo.opengeo.org/geoserver/wfs",
                    featureType: "tasmania_roads",
                    featureNS: "http://www.openplans.org/topp"
                })
            });

            map.addLayer(wms);
            map.addLayers([wms, layer]);
            map.setCenter(new OpenLayers.LonLat(146.7, -41.8), 6);
        }
    </script>
</head>
<body onload="init()">
    <div id="rcp1_map" style="width: 100%; height: 100%;">
    </div>
</body>
</html>


texts: Updating/replacing GeoTiff data store on file system causes GeoServer to halt

texts:
I have a few hundred raster layers that are updated at least 4 times a day.

The GeoTiff's are updated directly on disc.

I'm running it as a windows service and I can sell the wrapper app logging the
[Fatal Error] :1:1: Premature end of file.

The raster are in GeoTiff format, and it looks like a error to do with the files being updated causes GeoServer to fail and hang.

I cannot upload via the REST API as GeoTiff is not support,
I might be able get the wrapper to restart geoserver,

Also a smaller percentage of Shapefiles are also updated at the sametime, but I can tell from the logs it's more likely to do with the Raster layers

How can I  deploy raster data on windows, as a service, and to be able to update the raster data without having to restart the geoserver service wrapper?


texts: Splitting line every x meters using QGIS?

texts:
I have a line that I ultimately want to split up into points. The points should be every 100 meter along the line. So I don't want to extract the nodes.

Are there any Open Source (QGIS, Python) tools around for that?

The use case is that I have a bus line without bus stops. Though I know every 100 meter the bus stops. This way I want to generate bus stops to use as a GTFS feed.


texts: Converting personal geodatabase feature class to coverage using ArcObjects?

texts:
how to convert  personal geodatabase fatureclass to coverage using arcobjects with vba

pls provide me any documentation or samples.


texts: How to make graduated symbology icons visible in the exported map images of MapServer of ArcGis?

texts:
We created a new mxd file. We created a new layer using "Make a Feature Layer" command. Then we selected the new layer and got into the Symbology properties. There we created Quantities > Graduated Symbols for some database field.

The Symbology icons are visible in Arcmap:

We published the map as a MapServer.

The Symbology icons are visible in the published Legend of the Arcgis MapServer as seen here:

But the Symbology icons are not visible in the rendered map image of the MapServer:

Is there a special setting one needs to do in order to make the symbology icons visible in the exported map images of the MapServer?


texts: Arcpy: How to use da.cursors to 'pivot' one table into another

texts:
Let say I have a table that looks this:

State--City
TX--Houston
TX--Dallas
OK--Stillwater

And I want to insert it into a table that looks like this:

State--City1--City2--City3

For results that look like this:

State--City1--City2--City3
TX--Houston--Dallas--Null
OK--Stillwater--Null--Null

Should I use a search cursor to build a list then use that list to populate the 'insertRow'?


texts: Getting a KML file from embedded Google map

texts:
I need to get the coordinates for the markers on this map at the website :
https://www.evm.co.il/%D7%9E%D7%A4%D7%AA-%D7%A2%D7%9E%D7%93%D7%95%D7%AA-%D7%98%D7%A2%D7%99%D7%A0%D7%94-%D7%9C%D7%A8%D7%9B%D7%91-%D7%97%D7%A9%D7%9E%D7%9C%D7%99/
I already tried enlarging the map and getting the kml file but this map is embedded and not a my map format.
I also tried to get a json file by doing these steps :
open chrome developer (F12) > go to network tab > refresh the page > go to XHR > no json was found there.
I already checked with the website author and they don't mind if I use the map for my school project.
Is it even possible to extract the information from the website like this?


texts: Converting projected coordinates for Botswana

texts:
I live in Botswana.
How can I convert the coordinates they gave me on a lease?
I think it is very old system they use. Google Maps do not accept the coordinates. and ?earthplane? web age is only for USA.

On the lease is:

A=Y59665 X2493871
B=Y59460 X2493980
C=Y59701 X2494360
D=Y59819 X2494297

How do I convert the X & Y to degrees, minutes and seconds?

I have been searching in every Google site for a converter for these types of coordinates.


texts: Why am I losing part of my raster when I reproject?

texts:
QGIS 2.0
Working with tiff raster file.
Reprojections are in WGS84 Mercator.

Here is my problem: I am trying to use buffers or urban areas to analyze light intensity on a raster layer. I want 30 miles, 60 miles, and 100 miles for the buffer sizes. I have reprojected the vector layers to WGS84 / World Mercator. When I ran a zonal statistics, everything went smoothly until I was double-checking the output. I saw "null" on a few of the rows of information in the data table. When I looked at the raster image, which I also reprojected to WGS84 / World Mercator, a large chunk of the raster layer was missing. Images below.

Non-reprojected raster layer:

Reprojected raster layer:

I am assuming there is a problem with my reprojections. I am still new at GIS, so please bear with me. Thanks in advance.


texts: Handling errors from ArcGIS Feature Service applyEdits operation

texts:
When users calls applyEdits on our FeatureService, and get's something wrong with the request - it just returns "unable to complete operation".
Is it possible, to somehow get the detailed error?
I know that detailed errors are logged in ArcGIS logs.
I want to see those errors and return them to the client.
I can build SOI that could handle those requests, that go to the applyedits, that would be the preferable place to handle this. (We are using ArcGIS server 10.6.1)
Option that states "build your own SOE that does what applyedits does" is not acceptable, unless it is "easy" to build.
We had our own save operations made from times, when there were no applyedits. To increase the performance and usability, we are migrating our applications to use FeatureService applyedits.
Specific reason - why - we are having some DB triggers, that does some validation. And they throw detailed exceptions. (These work also, when user use ArcMap - so they are important).
Feature service just returns "Unable to complete operation"....
I also asked this question in Esri community forum.


texts: How do you delete temporary graph files using Arcpy?

texts:
I'm using the arcpy.MakeGraph_management and arcpy.SaveGraph_management  functions to create & save graphs in a Python Toolbox. The code works initially; however, subsequent use of the tool results in an error saying that the output from arcpy.MakeGraph_management (i.e., "temp_graph") already exists. If I close ArcMap and reopen it, I can again use this code fine, suggesting that the output from arcpy.MakeGraph_management is temporary; however I can't figure out how to delete this temporary graph in order to run the tool in succession or create multiple graphs in the same script.

# User provided input
folder = parameters[0].valueAsText
name = parameters[1].value
arcpy.env.workspace = folder

# Reference files
referenceFiles = os.path.join(os.path.dirname(__file__), "Reference Files")
graph_template = os.path.join(referenceFiles, "template_graph.grf")

# Create graph
arcpy.MakeGraph_management(graph_template, [...], "temp_graph")
arcpy.SaveGraph_management("temp_graph", "{0}/{1}_graph.pdf".format(folder, name))

Attempts at solution so far:

1)

arcpy.MakeGraph_management(graph_template, [...], "temp_graph")
arcpy.SaveGraph_management("temp_graph", "{0}/{1}_graph.pdf".format(folder, name))
arcpy.Delete_management("in_memory")

2)

if arcpy.Exists("temp_graph"):
    arcpy.Delete_management("temp_graph")
    arcpy.MakeGraph_management(graph_template, [...], "temp_graph")
    arcpy.SaveGraph_management("temp_graph", "{0}/{1}_graph.pdf".format(folder, name))
else:
    arcpy.MakeGraph_management(graph_template, [...], "temp_graph")
    arcpy.SaveGraph_management("temp_graph", "{0}/{1}_graph.pdf".format(folder, name))

Yet these tries still get me the same error:

I found someone with the same issue posted years ago: https://geonet.esri.com/thread/55685

I posted a response there in addition to posting it to Stack Exchange.

Edit (5/11/17):
I've discovered that the tool will run successfully if I set arcpy.env.overwriteOutput = True, although I'd prefer an actual solution with deleting "temp_graph". I had tried @Midavalo's suggestion to use arcpy.CreateUniqueName, which was unsuccessful, but forgot to include it in the original post.

I've noticed that despite my attempt to setup a Current Worskpace, it's not working.:

This is so despite the code where I attempt to set my current workspace and also try using arcpy.Describe() to print the path of where "temp_graph" is located.

param0 = arcpy.Parameter(
    displayName="Output Folder",
    name="in_workspace",
    datatype="DEWorkspace",
    parameterType="Required",
    direction="Input")
...
folder = parameters[0].valueAsText
arcpy.env.workspace = folder
arcpy.AddMessage("current workspace is {0}".format(arcpy.env.workspace))
desc = arcpy.Describe("temp_graph") OR desc = arcpy.Describe("{0}/temp_graph".format(folder))
print(desc.dataElement.catalogPath)


texts: Understanding garbage hillshade output from ArcGIS for Desktop?

texts:
I've used Hillshade successfully in the past, but for this DEM I keep getting a weird output. It's a 30-meter TIFF downloaded from the National Elevation Dataset.

And this is what the output looks like:

Kind of cool, but not at all correct. I've done other Surface analyses with the same DEM (e.g. Viewshed) and they ran without any obvious problem. Any ideas what may be causing the problem, and/or how to fix it?

Data came in UTM, so I re-projected it to the appropriate South Carolina State Plane. (I also copied it from the original TIFF into a geodatabase, ESRI raster.) It took a lot longer to process the same area, which got my hopes up, but the results are a similar "mesh" of gray squares/lines. I'm going to try hillshading another tile from the same source and see if it has the same problem.

The tile to the south hillshades just fine (without reprojecting or anything), and processes very quickly, so this seems to indicate there's something specific to the tile I want that is corrupted (or otherwise borked). It was downloaded from the NRCS Geospatial Gateway -- South Carolina, Richland County, 30m DEM, and it's tile ned30m33080.tif. (This FTP link should work for the next few days.)


texts: Converting encoding using ogr2ogr on macOS

texts:
I'm trying to convert my shapefiles from ISO-8859-10, also known as latin6 using ogr2ogr on my Mac like this (learned from here):

export SHAPE_ENCODING="latin6"
ogr2ogr file_utf8.shp file.shp -lco ENCODING=UTF-8

But, nothing is converted. Well, a new shapefile is created but it still has the same old encoding. I've tried with SHAPE_ENCODING="ISO-8859-10" and still no happy result. Then I tried just a random word like SHAPE_ENCODING="blablabla" but then ogr2ogr gives me an error:

Failed to create field name 'FIELD1' : cannot convert to blablabla

So it seems like ogr2ogr understands the latin6 encoding but then doesn't do the job for me. Maybe this is a macOS only problem?

If I do the conversion to GeoJSON with SHAPE_ENCODING="latin6" the GeoJSON file gets recoded to UTF-8 fine:

export SHAPE_ENCODING="latin6"
ogr2ogr -f GeoJSON file_utf8.geojson file.shp

Then I can convert the GeoJSON file back to Shape again:

export SHAPE_ENCODING="utf-8"
ogr2ogr -f "ESRI Shapefile" file_utf8.shp file_utf8.geojson

Is there a way of converting encoding without going through the GeoJSON conversion process?


texts: Finding distance between consecutive points using ArcGIS Desktop?

texts: I am wondering how I can get the distance between consecutive GPS locations in ArcMap 10 (i.e. the distance between consecutive GPS locations sent in during a wildlife telemetry study).


texts: How can I work with 3+ Gigabyte TIFF Files?

texts:
While I don't work in GIS myself I'm trying to support one of our users here who looks after all our mapping and GIS data. We've received a 3.2GB Tiff image containing satellite imagery of our site and I'm wanting to find a way to re-size it down a bit to make it easier to work with.

Is there any open source/free software out there that can handle this sort of file size?

I've tried Irfanview, Imagemagick, GIMP without success.

The servers I'm trying to work on have plenty of memory and disk available but perhaps not enough to fit the whole file.

Is there an easy way?

Thanks.


texts: I would like to have a logarithmic scale legend for my raster, anyone know how to do this in QGIS?


texts:
New to this, but played with Mapbox.

My interest lies in creating really good looking slippy maps, I am a web developer so I am interested in both the creation and interacting with the data using a API or similar.

So far I have created some lovely detailed building maps in Mapbox and been able to interact with the later data using MapBox GL - all good but I will soon surpass the limit and need to pay $499 a month my institution wont front that kind of money!

As Mapbox is open source and by definition available is it possible to re-create the project using Mapbox Studio open source and host on my own server?

I see the GitHub repository here https://github.com/spatialdev/mapbox-studio, but many of the links are 404 not found and it doesn't seem well supported - last commit was in 2015!!


texts: How to use d.* from a python/bash script in GRASS GIS 7.0?

texts:
Do d.* commands from GRASS GIS 7.0 run from a python script? I'm trying with this code:

import grass.script as grass

grass.run_command('d.mon', start="wx3", resolution = '1')
grass.run_command('d.mon', select="wx3")
grass.run_command('d.rast', map="M01relief")
grass.run_command('d.mon', select="wx3")
grass.run_command('d.vect', map="M03PTEtramos", color="black",  width='2')

When executing the complete script, the wx3 monitor starts, but remains empty. However, if I paste the commands one by one in GRASS Python Command-line, I get the layers drawn in the monitor. And using Cairo driver, I also get the layers drawn in a PNG file (directly from the script).

-

PS1: In the aforementioned Python code, the monitor seems not to apply the code, because if I check for the monitor commands, I get no info:

grass.run_command('d.mon', flags="c")
[Empty output]

(The monitor wx3 opens and remains empty)

Using a Bash script, the monitor is again empty:

d.mon start=wx3 resolution=1
d.rast map=M01relief
d.vect map=M03PTEtramos color=black width=2
d.redraw

But if I ask then for the list of commands I get results:

d.mon -c
Enlistar comandos para el monitor <wx3>:
d.rast map="M01relief" bgcolor="white"
d.vect map="M03PTEtramos" layer="1" display="shape" type="point,line,boundary,area,face" color="black" fill_color="200:200:200" width=2 width_scale=1 icon="basic/x" size=5 label_layer="1" label_color="red" label_bgcolor="none" label_bcolor="none" label_size=8 xref="left" yref="center"

Actually, doing d.redraw from command line (not from the script), shows the map.

The problem happens both in Linux Mint 64 bits and Windows7 64 bits (the second running in a VirtualBox machine)

PS2: This is weird. The commands are written to the cmd file of the monitor (e.g.: .tmp/24518.0.cmd):

d.rast map="M01relief" bgcolor="white"
d.vect map="M03PTEtramos" layer="1" display="shape" type="point,line,boundary,area,face" color="black" fill_color="200:200:200" width=2 width_scale=1 icon="basic/x" size=5 label_layer="1" label_color="red" label_bgcolor="none" label_bcolor="none" label_size=8 xref="left" yref="center"

But I do not get the layers shown until I manually open the file (in a text editor) and save it. Using some python code to open, edit and save the trick doesn't work [here 'rutamon' is the path to cmd file of wx3 monitor]:

with open(rutamon, 'a') as comandos:
    comandos.write('d.redraw')
    comandos.flush()
    os.fsync(comandos.fileno())

PROVISIONAL WORKAROUND: I have found that if I open a second monitor after the first, the layers are shown:

grass.run_command('d.mon', start="wx3", resolution = '2')
grass.run_command('d.rast', map="M01relief",  overwrite="true")
grass.run_command('d.mon', start="wx0", resolution = '1')
grass.run_command('d.rast', map="M01relief",  overwrite="true")
grass.run_command('d.mon', select="wx3")
grass.run_command('d.redraw', verbose="true")

With this code, I get an empty "wx0" monitor, and a "wx3" monitor showing my layers. But I need that second 'ghost' monitor, which is not very elegant... Any help with this is welcome.


texts: Finding GIS data sources for developing nations?

texts:
Not all countries have the benefit of funding/technologies for developing GIS data stores. I'm looking for a list of common GIS data sources that should be investigated when trying to build a good basis for GIS analysis in developing nations (e.g road/river networks, DEMs, land use). Free or paid data.

My first stop is usually Glovis for Landsat coverage, but I often come up short on vector data.

For example, I would like to create a series of thematic maps for Nepal containing any and all GIS information I can get my hands on. Where would you begin?

Organizational sources could also be useful, even if their data is not available online. E.g contacting the government department in charge of water resources.


texts: Is it possible to change orientation of specific pages in QGIS print composer?

texts:
I noticed I can add multiple pages to a print composer - cool!

Can I use an expression to set the orientation of specific pages?

(I know a little python so if a function is more appropriate I'd be open to exploring that!)

I have about 5 pages, but need pages 2 and 4 to be landscape, whereas pages 1,3 and 5 should be portrait.

There is a layout_numpages variable, but that's only for the number of pages, but in the Item properties tab, there is a layout_page variable under Composer item.

So I tried forcing the orientation with an expression like:

CASE
  WHEN  @layout_page IN (2,4)
    THEN 'landscape'
    ELSE 'portrait'
END

But that didn't have any effect. The other thing I can do is rotate all the map elements 90deg, but that's a little more work than just rotating the actual page to match the map elements.


texts: Updating print composer map when template in new project is used?

texts:
I have a print composer template saved, but when I open it in a new project I can't see a way to get the map to update to the new project. It just stays blank.
Is there something like the legend where you can auto-update it or link it to the current project or do I need to add the map in manually every time?


texts: How can I change object image in MapInfo using MapBasic?

texts: I need to change the image of the objects which are in the Selection Table using MapBasic functions. When I say "Image" I mean point object icon on the map.


texts: Computing Layer statistics using PyQGIS script

texts:
I am trying to develop a script for the Processing Toolbox which computes layer statistics, creates a style of upwards pointing triangles for higher than average elevation and downwards pointing red triangles for lower elevation based on statistics of a vector layer. This is how far I have gone:

##GeoGang=group
##GeoGeeK=name
##Input_vector_layer=vector
##Keep_attributes=boolean True

The output


texts:
So I am trying to batch save 100's of styled rasters. I have managed to save a single raster in the python console using the following code:

layer = iface.activeLayer()

extent = layer.extent()
width, height = layer.width(), layer.height()
renderer = layer.renderer()
provider=layer.dataProvider()
crs = layer.crs().toWkt()

pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())

file_writer = QgsRasterFileWriter('e:/Data/output.tif')

file_writer.writeRaster(pipe,
                        width,
                        height,
                        extent,
                        layer.crs())

What do I need to add so I can save all the rasters currently open in QGIS?


texts: Filling in new field in attribute table based on values of another field in ArcGIS Desktop?

texts:
I have added a new field (Category) in an attribute table of a shapefile and I would like to fill it with numbers, 0 or 1, based on the numbers of another field (Calorific value) of the same attribute table.

That means that if the Calorific value is higher than 15, then I would like to fill in the corresponding cells of the (new) Category field with 1. If the Calorific value is less than 15, then the rest cells of the (new) Category field to be filled in with the number 0.

Is there a way to do this automatically rather than manually?


texts: ArcObjects Add-in - Making distinct add-in components talk to each other

texts:
I'm starting developing with ArcObjects and I'm stuck at trying to make two comboboxes talk to each other. What am I trying to do is changing the contents/selection of ComboBox B based on the selection of ComboBox A. But so far I don't know how to get the ComboBox B instantiation inside ComboBox A code.

An example code I found so far use ICommandBars, but it needs an IApplication reference that ComboBoxes do not receive as parameter in their initialize() method. Another example is by using method findExtensionByCLSID(UID), but it also needs an IApplication reference.

Can someone give me a hint on how to grab another add-in component informations (such as selected item in ComboBoxes)?

A piece of code I have to try catching a ComboBox instantiation is presented below (this code is inside a combo named DimCombo, and tries to pick a FormaCombo combo instance to do some work). I guess my problem is at setting the app variable, as it is not given in combos initialization method.

public void getFormaCombo(){
    try{
        UID pUID = new UID();
        pUID.setValue("br.mil.eb.dsg.arctools.formacombo");
        formaCombo = (FormaCombo) app.findExtensionByCLSID(pUID);
    }
    catch (Exception e){
        e.printStackTrace();
    }
}

Note: br.mil.eb.dsgtools.formacombo is the ID of the desired combo instance (set in config.xml). I'm using Java, but hints can be in any language.

Edit

Apparently I found a way to grab ComboBoxes from within a Tool code. Solution is a code similar to the previous one, but it is now included inside the desired ComboBox implementation. Client tool should invoke this method to fetch the desired combo's instance.

public static DimCombo getExtension() throws UnknownHostException, IOException{
    if(s_extension == null){
        UID pUID = new UID();
        pUID.setValue("br.mil.eb.dsg.arctools.dimcombo");
        app.findExtensionByCLSID(pUID);
    }
    return s_extension;
}

Tool components need to do something like:

DimCombo combo = DimCombo.getExtension();

And then ComboBox methods are made available (like getSelected()).

But I am still with no solutions on making another ComboBox call the getExtension() like above. It gives a null pointer exception exactly at the above line.

So I still need some more help.


texts: Adding polygons to existing attribute table where data is stored as points?

texts: I have a shapefile of golf courses saved as points, I would like to digitise the boundaries and add the polygons to the existing attribute table.  I have opened editor and clicked on the point data shapefile but have no option to create a polygon.


texts: Attribute Table fields NULL for volume calculation - confused

texts: I just tried a volume calculation with QGIS on a DSM. After using the zone statistics function the attributes (mean, max, min) for the corresponding polygon are all NULL. Anyone knows why that doesn´t work?


texts: How to do a POST instead of GET when calling sql.execute()?

texts:
I am calling the CartoDB sql.execute() function with a rather large WKT geometry in the SQL query and this is causing a "414 Request-URI Too Large" error on the server b/c it is making a GET request.

My question is how to I make the sql.execute() function make a POST instead of a GET so I can perform a large SQL query?


texts: GeoServer - layer is read only (it already has a primary key)

texts:
I created a simple table in my Postgres\PostGIS DB. I then created the GeoServer layer with the SQL:

select * from tablename

I'm trying a WFS Transaction, like I usually do (and it works on similar tables!) to add data to that layer. I get this error:

{http://myorganization.org}layername is read-only

It's a common error, that's usually due to the lack of a primary key in the "tablename" table. But the table has a primary key!

create table tablename(
serialID serial PRIMARY KEY,
field1 text,
field2 text,
field3 text,
field4 timestamp without time zone,
geom geometry
);

I already have other tables where the Transaction works in the same GeoServer instance, so it's not a system wide option or lack of permission to the tomcat user. How do I fix this error?


texts:
I am working on a data collection application for Android and I am able to use CartoSQLService to read information from my table in Carto, however I am having difficultly adding data. I have searched online for any guidance with no luck.

Does anyone know the proper procedure to add data to a Carto table using the Carto Mobile SDK? Is the API Key suppose to be included? If so, how is this done since the Carto documentation does not show how this is accomplished for Android applications?

Here is what I have so far:

final CartoSQLService service = new CartoSQLService();
service.setAPITemplate("https://{username}.carto.com");

Thread thread = new Thread(new Runnable() {
@Override
     public void run() {

          try {
            service.queryData("INSERT INTO {table_name} (lat,long,sod) VALUES (" + locArray[0] + "," + locArray[1] + "," + checkBox.isChecked() + ")");

        } catch (IOException e) {
            e.printStackTrace();
        }
    }});
thread.start();

Error:

E/carto-mobile-sdk: HTTPClient::makeRequest: Bad status code: 401, URL:

{
    "error": [
        "permission denied for relation acorn_table"
    ]
}


texts: What are the top 5 file types for 3D model?

texts:
I don't know much about GIS but I'm investigating a project that requires a 3D model of a hilly large town (town centre, suburbs, golfcourse, etc.).  I'm not certain which file-type of these Geospatial Formats will give me the best data for the model.

I need the data to be able to distinguish between homes and their lawns (roof height), spaces between town buildings (most are less then 10 stories; with roof height), hills, valleys, and so on.  Perhaps a resolution of 0.5 cubed-meter or less; I am guessing.

I need to create a 3d model of the surface of the area where the surface includes the buildings, trees and the like.  I'll then be doing a few flow simulations on this 3D model.  The is seeing how water (flood, industrial run-off), industrial smoke stack soot (pollution), industrial noise, and the like, travels through the area.

I plan on using the Mathematica technical computing software (+ some cold computing resources) for this and it only has the few Geospatial Formats mentioned in the link. It can make 3d models from point clouds and from elevation maps.  I'm most concerned about the resolution I'll be able to get as I need the model to be able to make out buildings from other surrounding buildings in both in the recreational, residential, commercial, and industrial areas. Large tree line heights would also be very good to have in the data.

The data will be from an aerial survey but I need to specify the file-type.


texts:
I would like to make a map displaying how several building uses have changed during the last decades. For example, this building was built in the year 1483 as a monastery, in 1880 it turned into a seminary and since then it has been a prison (1936-1939), has been abandoned (1939-1940), then transformed into a boarding school (1940-1998) and currently is a camping house (1998-present time).

I was thinking of adding two new date fields start_date and end_date in the layer with the geometry and other attributes and then combine everything with TimeManager plugin so I could dynamically see the evolution. Unfortunately that's not enough, as 1) the uses may change whereas the geometry will be the same one; and 2) that would only allow me to have a single use per row instead of having a complete uses' history. That made me think that I need to have more than one row for the same building if I want to have a "history" of uses. For that reason I was then thinking of creating a layer for the geometry and a different one for the building uses with the following fields:

id
geometry_id
description
start_date
end_date

By doing this I could easily "relate" a single geometry with one or more uses (each of them being a different row, each one having its own start and end dates). Unfortunately I do not know how can I relate both layers in an easy and fast way (preferably using GUI - I don't think that doing a regular join would be enough, as I have no easy way to know the geometry_id values beforehand) nor I am sure that is the best approach to do so (looks somewhat complicated to manage if I have many buildings' that I want to keep track of their changes, as it is the case).

What would you recommend me to do?


texts:
I was wondering if there is a way of passing a variable as the XML query in mapnik.filter().

Namely, we know we can do the following:

rule=mapnik.Rule()
rule.filter=mapnik.Filter("[NAME]='test'")

to obtain a map which does something special if the NAME of something is 'test'. However, I would like to do the following:

s="[NAME]='test'"
rule=mapnik.Rule()
rule.filter=mapnik.Filter(s)

This, however, doesn't seem to be working, giving the following error:

Boost.Python.ArgumentError: Python argument types in
mapnik._mapnik.Expression()
did not match C++ signature:
Expression(class std::basic_string<char,struct std::char_traits<char>,class
std::allocator<char> > expr)

Here the OP works directly with Boost Python and is able to convert his string to an immutable one; I, however, am just working with Mapnik, which is a package which uses BoostPython in the background, so it's unlikely a similar solution would work for me.


texts: Trouble with PostGIS and multiple selection sources

texts:
I'm attempting to take geometry columns from two different tables, perform basic PostGIS operations on them, and have them returned in a single result set.

My SQL statement looks like the following:

SELECT (ST_Dump(ST_Polygonize(exterior_rings.the_geom))).geom AS exterior,
       (ST_Dump(ST_Polygonize(interior_rings.the_geom))).geom AS interior

FROM (SELECT the_geom FROM geodata."_unioned_exterior_rings") AS exterior_rings,
     (SELECT the_geom FROM geodata."_unioned_interior_rings") AS interior_rings;

When I execute this command in pgAdmin III's arbitrary SQL query tool, after roughly 8 minutes of execution, I receive:

********** Error **********

... and the query tool complains it is no longer connected to the database.  This leads me to believe the operation is timing out.   The respective sizes of the tables are fairly humble.  _unioned_exterior_rings has ~4000 rows, and _unioned_interior_rings has ~100.  I've had other queries succeed after twice as much run time so I don't feel as if I'm hitting some preset limit.

When I use pgAdmin's console tool to execute the query, after some time it becomes apparent that the entire tool has been disconnected without any kind of error message.  On some reconnect attempts I've received notices that the reconnection fails because the database is in a fatal recover mode.  If I try again shortly there after it usually works.

When I execute the statements separately (doing the polygonize/dump on an individual table) I get the expected results after only a few seconds.

The version of PostgreSQL is 8.4, PostGIS is 1.4.2, and pgAdmin III is at 1.12.2.

I can't claim any real SQL proficiency, so I'm willing accept that my mistake may be simple or fundamental.  Any help would be appreciated.


texts: Failed DEM via Natural Neighbor in QGIS 2.14.2-Essen

texts: Following the indications of Monde Geospatial (https://www.youtube.com/watch?v=bLbY3iMBW-A) I created a CSV file with the altitudinal data (Google Earth) for 2182 total points and imported it into QGIS 2.14.2-Essen. I am stalled at the next step: DEM generation. The message "The following layers were not correctly generated. Grid" appears on the Natural Neighbour log console. Before the failed output and as the process is underway, I observed a "Can´t load requested DLL" appearing in the same console. I need help identifying the cause and, if necessary, seeking alternatives to Monde Geospatial´s method.


texts:
So I have second level administration data of a country with shapes in a PostGIS database. Now I actually need the first level administration level representation of the country, but here I only have the name of that first level admin and the names of its subdivisions as a text file (which do not need to match completely with the names I have in the DB).
Now I would like to manually remove the admin2 data entries and replace them with the admin 1 data. However, for the geometry column I was thinking to do a ST_UNION by ids of the admin 2 areas of one admin1 area. Is this possible? Right now there is no indicator of shared admin 1 information in the DB. Something like:
ST_UNION((SELECT geom FROM states WHERE id IN (1,2,3,4)))


texts: Differential correction wizard doesn't pop up?

texts:
Any Trimble GPS Pathfinder Office users out there ever run into this problem? When I try to run the differential correction wizard the program stutters for a brief second then nothing happens. All other utilities (data transfer, export, etc.) still work fine.

Also, if I decide I need to reinstall the program does anyone know if I have to deactivate my single use license beforehand?


texts: gdal2tiles_parallel seems to be signiciantly slowing down for super high-res images

texts:
I have 3 images that I've been trying to 'tileize'

1 86,000*41,000
2 86,000*21,000

both are simple .pngs

I've been using the parallel version of gdal2tiles on a DigitalOcean box with 20cores and 64GB of Ram(Ubuntu 14.04 x64) with the package gdal-pin installed.

Whenever I try and tileize the images with the command

python gdal2tiles_parallel.py -v -p raster -k Something.vrt -z 1-7 -r lanczos

it starts off fine but then after a minute or two it just grinds to slower than dial-up snail speed. A few images a minute. Is there a configuration option I'm missing or is the image simply too big? The ones I'm going to be processing are 800mb to 1.5gb.


texts: Cliping a Raster to a polygon producing errors

texts:
I'm clipping a raster to a polygon using extract by mask (Arc 10.5). This produces a very skewed result. White washed image. I've tried using raster clip, but the "use input features for clipping geometry" option is grayed out so that I can't check it meaning I'm not getting the result I want. Some details: when I enable DRA the image looks much better. But I need to be able to view it properly without clicking enable DRA.

Any insight as to why the clipping geometry option is grayed out would be great. I thought it might be a memory issue but I have deleted a lot of data and still no change. Prior to the clip the image is a landsat 8 16 bit image with compression type LZ77. After the extract by mask the image becomes compression type LZW. The image that worked properly was compression type none, but when I export the data and select compression type none the image somehow becomes 32 bit.


texts: What does percentile mean in simpleComposite of Google Earth Engine?

texts: I was using 'ee.Algorithms.Landsat.simpleComposite' to composite images, but I'm not sure what does this 'percentile' parameter mean?


texts: Input string value in QGIS 3 Processing Scripts

texts:
I am trying to port my old QGIS 2.x processing scripts to QGIS 3.x.  One of my scripts just needs two variables to initiate several processes to set-up a working environment in my QGIS map document (downloading layers, subsetting data, etc).  Using the template script (QGIS 3.2.0) I have changed the static variables, 'INPUT' and 'OUTPUT' to 'START' and 'END'.  I then changed the initAlogrithm function to:

    def initAlgorithm(self, config=None):
    """
    Here we define the inputs and output of the algorithm, along
    with some other properties.
    """

    self.addParameter(
        QgsProcessingParameterString(
            self.START,
            self.tr('Model Start Time'),
            'YYYYMMDDHH'
        )
    )

    self.addParameter(
        QgsProcessingParameterString(
            self.END,
            self.tr('Model End Time'),
            'YYYYMMDDHH'
        )
    )

Then in the processAlgorithm() function, I tried calling the inputs as:

start_time = self.parameterAsString(
        parameters,
        self.START,
        context)

This yields the following error:

TypeError: invalid result from AssimStartUp.processAlgorithm(), NoneType cannot be converted to a C/C++ QMap<QString,QVariant> in this context

So, what am I doing wrong that is prohibiting the string from the input window being passed into the processAlgorithm() function?


texts: Is editing geometric network in versioned geodatabase any different from when such data is stored elsewhere?

texts: I need to work on "Checked Out" version of a geodatabase, editing geometric networks. What will be the difference in terms of the working process of the editor ?-do i need to use some special toolbars ? is it possible to work on versioned GDB with ArcEditor 10.1 ?


texts: OpenLayers 2 - Capture failure of WFS HttpRequest

texts:
I have a WFS Layer whose HTTP Request might fail from time to time. How do I capture the failure?
        var opt =  {
            url: widget.options.wfs.url,
            params: {
                request: "GetFeature",
                service: "wfs",
                version: "1.1.0",
                typeNames: widget.options.wfs.featureType.split(",").map(function(str) { return str.trim(); }),
                srsName: widget.map.getProjection(),
            },
            format: new OpenLayers.Format.GML({})
        }

        var wfsLayer = new OpenLayers.Layer.Vector("WFS", {
            strategies: [new OpenLayers.Strategy.BBOX()],
            protocol: new OpenLayers.Protocol.HTTP(opt),
            projection: widget.map.getProjection(),
            styleMap: getStyleMap(),
        });

        widget.map.addLayer(wfsLayer);


texts: Displaying icon callout when no text in data field using QGIS

texts:
I am using QGIS.
On my maps, I have SVG Parking symbols. I need to use callouts to move them away. I have also a reference number beside many of these icons, but not all.
How best to set this up?
I tried this and it almost works...
In the layer styling background tab, I added my SVG symbol. Created a field column for my reference #s and it displays. The original location marker point I do not need to see, so I disable the symbol layer.
What I cannot solve is how to display the icon when there is no reference # in the field table, just NULL?
Since it needs to see text in that field to display the background parking symbol, how do I find a workaround?
Is there a trick or is my approach wrong?

To elaborate more on this issue, here are some screen caps. My Parking Icon with the reference number beside it as a callout. I like this no problems.

Here are my setting in the Background tab. I display the blue square P marker using the Background tab becouse I do not know how else to get the callout feature to work.

In the main Symbology window I turn off the layer that would show a marker at the original point. I do not what this to show as I will be scaling the maps and all those markers would colide on my rendered maps.

Lastly here is my Parking data table with my PLot column that I added. Note how fid 156 displays P4 because there is the number 4 in the field. But row 293 does not display P becouse there is no number in the field. (And I do not want one associated with that marker.)
Yet I do what it to show up just as a P marker with no number beside it. How?


texts:
I want to create a vrt from my plugin via gdal-tool. From OSGEO4W-Shell the command works. In my code i tried this:

QString program = "gdalbuildvrt";
                    QStringList arguments;
                    arguments.append("C:\\Temp\\output.vrt");
                    arguments.append("P:\\Daten\\Rasterdatenkatalog\\Raster1\\kat_1\\teil1\\4331ein.tif");
                    arguments.append("P:\\Daten\\Rasterdatenkatalog\\Raster1\\kat_1\\teil1\\4332ein.tif");
                    arguments.append("P:\\Daten\\Rasterdatenkatalog\\Raster1\\kat_1\\teil1\\4333ein.tif");

                    QProcess p;
                    p.start(program, arguments);
                    if (p.waitForStarted())
                    {
                        p.waitForFinished();
                        qDebug() << p.readAllStandardOutput();
                    }
                    else
                        qDebug() << "Failed to start";

I get no errormessage but the file will not be created. Do I need the full path to gdalbuildvrt? I thought in qgis the path is even set.

Thanks for any idea.


texts: Dynamic tools for ArcGIS python toolbox: Class factory works in python, but ArcGIS Pro only shows last tool/class created

texts:
I started with the the standard python toolbox template for ArcGIS: https://pro.arcgis.com/en/pro-app/latest/arcpy/geoprocessing_and_python/a-template-for-python-toolboxes.htm
My goal was to find a way to create tools/classes from a list of tool names dynamically using a base class/tool as a template.   Basically, the only difference in these tools is the name (label) and description text.  The parameters, execution and all other functions besides init should be shared/copied into each tool.
I have tried deepcopy and class factory methods for inheritance, but nothing has worked within ArcGIS Pro.  In python, the classes/tools seem to be made as expected and have the desired properties.  However, in the ArcGIS Pro toolbox, only the last class in the list displays.  I will give the example code below.   My main question is, should this work as is with ArcGIS Pro or am I missing something with the return classes.
import arcpy

class baseTool(object):
    def getParameterInfo(self):
        #Define parameter definitions
        p0 = arcpy.Parameter(
            displayName='Some Value Input',
            name='SomeVal',
            datatype='GPBoolean',
            direction="Input",
            parameterType='Optional',
            enabled=True)
        p10value = False

    params = [p0]
    return params

def isLicensed(self):
    #Set whether tool is licensed to execute.
    return True

def updateParameters(self, parameters):
    #Modify the values and properties of parameters before internal
    #validation is performed.  This method is called whenever a parameter
    #has been changed.
    return

def updateMessages(self, parameters):
    #Modify the messages created by internal validation for each tool
    #parameter.  This method is called after internal validation.
    return

def execute(self, parameters, messages):
    #The source code of the tool.
    # do some standard process
    arcpy.AddWarning(self.label)
    arcpy.AddWarning(self.description)
    return

# class factory method

def classFactory(name, desc, BaseClass=baseTool):
    def __init__(self):
        #Define the tool (tool name is the name of the class).
        self.label = name
        self.description = desc
        self.canRunInBackground = False

    newTool = type(name, (BaseClass,),{"__init__": __init__})
    return(newTool)

# list of tools to create from template/baseTool
inputTools = {  'Tool1': 'Some text info',
                            'Tool2': 'Some text info',
                            'Tool3': 'Some text info',
                            'Tool4': 'Some text info' }

# iterate the list to create a new tool for each
toolObjs = []
for k, v  in inputTools.items():
        newTool = classFactory(k, v)
        toolObjs.append(newTool)

class Toolbox(object):
    def __init__(self):
        #Define the toolbox (the name of the toolbox is the name of the.pyt file).
        self.label = "Some Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = toolObjs


texts:
Is there any way to capture ArcMap Drawing Errors inside of a python script?

Through python, I would like to iterate through each layer inside of an MXD, make the layer visible, refresh the table of contents and the active view. I then want to see if any of these layers cause a drawing error, in which case I can write the layer that caused the drawing error to a log file so that I can later review and fix the problem causing the drawing errors. The only problem is I have no idea of how to capture the message from the ArcMap Drawing Error Dialog with python. I don't even particularly care what the error message is, I just want to see if an error is generated.

Example of an ArcMap Drawing Errors Dialog:


texts:
How to use symbolizer in openlyers3? I need to display my line with double dashed.

We have option in 2 like symbolizer : {strokeDashstyle:"dash"}

Please help me..


texts: Shortcut to display scale from cache in qgis print composter

texts:
When loading a map in the print composer I always click on the main area where the map displays, then click on 'set to map canvas extent'. Then I adjust the scale as I need. This is all fine.

How do I pick up the scale (e.g. 4000) and have this value automatically inserted into a different part of the printing template as follows?

  Scale 1:4000 @ A3

At the moment I have to type the scale by hand every time.


texts:
I am trying to export a DEM raster into a xyz file (for later use).
I first clipped the raster in order to limit its extent and the lines in the xyz file.
However, gdal2xyz does not ignore the nodata values (set to 0 here but it is the same with -9999) resulting in a very large .csv files, with most lines with a Z of 0 (useless in my case).
Is it possible and how to ignore those nodata values with gdal2xyz in order to reduce significantly the size of the file?


texts: QgsVectorDataProvider addFeatures doesn't take attributes

texts:
I'm trying to add features in empty QgsVectorLayer with addFeatures, I got the layers without attribute:

def create_features(input_file, wid, hite, focal_lgth):
    features = []
    filename_def = {
        'name': 'filename',
        'type': QVariant.String,
        'length': 50,
        'precision': 0
    }

    field = QgsField()
    field.setName(filename_def['name'])
    field.setType(filename_def['type'])
    field.setPrecision(filename_def['precision'])
    field.setLength(filename_def['length'])
    fields = QgsFields()
    fields.append(field)

    with open(input_file) as file:
        has_header = csv.Sniffer().has_header(file.read(1024))
        file.seek(0)  # Rewind.
        reader = csv.reader(file)
        if has_header:
            next(reader)  # Skip header row.
        for i, row in enumerate(reader):
            feat = QgsFeature(fields, i)
            ----------------------------
            feat.setAttribute("filename", row[0])
            feat.setGeometry(geometry)
            features.append(feat)
    return features

fts = create_features(input, 4000, 3000, 4.3)
layer = QgsVectorLayer("Polygon?crs=EPSG:4326", "temporary", "memory")
layer.dataProvider().addFeatures(fts)

While I printed fts[i].attribute("filename") I got the right value.


texts: Point features disappear when adding new data

texts:
Recently while adding features to a point layer i noticed that once i had saved my edits the already existing points disappeared. Although not visable on the map all of them and their data were still in the attribute table. Any guesses on how to make all of the points viewable again?

Thanks,

L


texts: Spatial join 2 polylines without common attribute

texts:
I need to spatially join 2 polylines, but the problem is that the tables have nothing in common. Also, one of the tables is 4 times as big as the other one. Essentially, what is need is to combine the 2 layers into one based on location. I tried with buffers, and simple joining didnt work, so I will try spatial joining, but im wondering is there is a surefire way to do the following:

A poorly digitized version (has what i need in the take)
Properly digitized version ( 4x as big as previou)

If a polyline in  #1 is found to be close to the polyline of #2, and share approximate angles, then it should merge.


texts: GDAL Python bindings not picking up OGR PostgreSQL driver

texts:
I have a windows server (Microsoft Windows Server 2003) with:

Python 2.7,
PostgreSQL 9.04-1,
GDAL 1.8.1 (installed from the windows installer
(gdal-18-1500-core.msi) hosted at
http://vbkto.dyndns.org/sdk/PackageList.aspx?file=release-1500-gdal-1-8-mapserver-6-0.zip), and
GDAL Python bindings installed from the windows installer
(GDAL-1.8.1.win32-py2.7.msi) at the same url as above.

If I do

ogrinfo PG:"dbname=postgis_tpl user=postgres password=xxxxxx"

connection with the database is successful:

INFO: Open of `PG:dbname=postgis_tpl user=postgres password=xxxxxx'
  using driver `PostgreSQL' successful.
1: world_merc (Multi Polygon)

However, from the Python shell, if I do:

>>> from osgeo import ogr
>>> print ogr.Open("PG:dbname=postgis_tpl user=postgres password=xxxxxx")

I get:

None

and:

>>> print ogr.GetDriverByName("postgresql")

gives me:

None

and if I list all the drivers available via the Python bindings:

>>> for i in range(0,ogr.GetDriverCount()):
...     print ogr.GetDriver(i).name

the following drivers are listed (PostgreSQL is missing):

ESRI Shapefile, MapInfo File, UK .NTF, SDTS, TIGER, S57, DGN, VRT, REC, Memory, BNA, CSV, NAS, GML, GPX, KML, GeoJSON, GMT, SQLite, ODBC, PGeo, MSSQLSpatial, PCIDSK, XPlane, AVCBin, AVCE00, DXF, Geoconcept, GeoRSS, GPSTrackMaker, VFK, PGDump, GPSBabel, SUA, OpenAir, PDS, WFS, HTF, AeronavFAA

My environment variables are (I think) set correctly:

GDAL_DRIVER_PATH = C:\Program Files\GDAL\gdalplugins

(There's nothing that looks like the PostgreSQL driver in the gdal-plugins folder anyway).

GDAL_DATA = C:\Program Files\GDAL\gdal-data

and

C:\Program Files\GDAL

is in my PATH (GDAL is installed at C:\Program Files\GDAL).

So what am I missing? Why can't Python see the PostgreSQL driver? Do I need to install its bindings separately?