texts: List layers from current view using Arcpy

texts:
I am using ArcGIS 10.4.1 and Python 2.7.10.

I have a mxd file of a rather large area which contain hundred of layers. I wish to work in a specific sub-area in that workspace, where only a few layers will be relevant.

I was looking for a way to output a list of the layers which are spatially overlapping with the current view (sub-area) in ArcMap. The layers could be visible or not, as checked (or not) in the workspace. The mxd is up to date with this view.

Is such a thing possible?

I am using as a basis the following code, which output all the layers in the workspace:

import arcpy

mxd = arcpy.mapping.MapDocument(r"C:\path\to\file.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
print arcpy.mapping.ListLayers(mxd,"", df)

EDIT: a practical example to illustrate would resemble to this:
Let's say I have dozen of layers from topographical maps and ground features (water, river, topoline-level, roads, electric line, settlements, glaciers, reservations, highways, airports, and forests) at a country scale.

I zoom in a very remote forested area, where only the layers forest, river, topoline-level are overlapping the view. All other are not relevant. I am looking if it is possible to arcpy those forest, river and topoline layers in a list and discard the other.


texts:
I have a MSSQQL Spatiatal layer with Point data (borehole locations, pointID) and a joint table lithology (soil info per borehole, pointID, depth from, depth to, constituent)

Is there a way to use the html map tip display to output a small table like this?

Depth from/      depth to/     constituent/
0/               1/            sand
1/               2/            clay
2/               3/            gravel

I manage to get the 3 columns but get only one row showing the last existing Lithology row for the chosen pointID ... goal is to see all existing rows for a specific pointID.

Any ideas??


texts:
I have the following Model, that runs the Location-Allocation of Network Analyst.

I'm having a problem, in the "Calculate Field Tool". First I add a new field (it's running good) and then I'm my model breaks, when it runs the "Calculate Field Tool".

The expression I'm using is:

New_field = %Name%

What I'm intending to do is to set the content of my variable "Name" retrieved from the Iterate Feature and write it in my new field.

However, I'm having a problem with that, resulting in a dialog box (see image below)

I checked the content of the "Name" variable and it contains the name I want and it's a text string. I don't know why it is not writing it to my new field (which is also a TEXT field)

Do you have a hint of what's going on?


texts: Using SRTM Global DEM for Aspect calculation?

texts:
I need to derive a raster representing aspect for North America (from STRM data, originally in geographic lat/long). In most projections, true north is only map north along eg the central meridian so aspect values are distorted.
I attempted calculating aspect in Mercator and then reprojectiong to eg Lambert Conformal Conic, but the aspects values are very pixelated and distorted.

The following answer by whuber suggests how slope could be calculated accurately across a continent: "Compute the slope (as rise:run or a percent) using a Mercator projection, then divide the result by the cosine of the latitude". Using SRTM Global DEM for Slope calculation?

How would this be extended to calculate aspect?

Would one use flow direction based on the corrected slope?

Or can other methods be used?

I see in Compute global DEM topographic indices and compromise projections that @whuber states "Having done this [corrected slope calculation], perform the followup calculation of flow directions, aspects, etc. using the corrected slopes."

I'm just not sure how to derive aspect from a slope raster (vs from a DEM) in ArcGIS or otherwise.


texts: Cancel event in Python Gp tools with arcpy

texts:
In GP-tools written in C# (ArcObjects) we can listen for the Cancel-event and react on that. I wonder if there is anything similar in Python GP-tools applying arcpy.
I googled that issue and the only thing I found was the Cancel-method on the Result-Object (so that you can ask the state of the result of a GP-tool, i.e. if it was canceled).
I'm aware that you can cancel GP-tools running in the background in the Result Window, but this Cancel-event is catched by the system and therefore the developer or the GP-tool cannot react on that.


texts: QGIS 2.0 How can i get a precisely Rendering time?

texts:
I work with QGIS 2.0.1.

I have a large number of raster data files and i selected different methods to represent them.

Method: Import all raster data as individual files
Method: Generate a single file with the help of the rastertool "Merge"
Method: Generate a virtual raster (catalog)

I would like to know if it is possible to get the exact Rendering time? For example how long qgis needs to render the map if you zoom in.

This would be very helpful to compare the rendering time between the three methods.
I know which method is faster but a time value would be important for me.

Thanks for your help


texts: How to merge rasters with the same but unequal cell size in ArcMap?

texts: I'm trying to merge raster tiles that all have the same, but unequal cell size (X: 2,8, Y: 2,9). When I use the Mosaic To New Raster tool, I can enter ONE cell size which either results in a raster shifted in X or Y-direction or a black line going through the resulting raster when I don't specify anything. Is there any way to get around this?


texts:
I am working on a script tool that will be run from a custom toolbox in ArcCatalog.

I have written everything with the inputs hardcoded into it for testing. Works great running from IDLE.

I had hoped that all I needed to do in the end was to remove the hardcoded paths and replace them for the input tool parameters:

workspace = arcpy.GetParameterAsText(0)
...
save_flag = arcpy.GetParameterAsText(10)

But when running from the script tool GUI, I'm finding several places that it fails.

I have several places where I use raster math in place raster calculator as shown below , but it errors as "not licenced" (I am sure have SA activated earlier in the script) or as "invalid parameters".

First attempt (works in IDLE, fails in script tool)

eopslopes = arcpy.Raster(masked_dem3) - (float(SinkVal) / (arcpy.Raster(euc_eop) + 1))

and shown below in context

arcpy.FeatureToRaster_conversion(eop, "OBJECTID", r"in_memory\eop_ras", dem)
inSourceData = r"in_memory\eop_ras"
slopeDist = 15
SinkVal = 1
buffdist = 10
cellSize = dem
arcpy.gp.EucDistance_sa(inSourceData, r"in_memory\eucdist_eop", "", cellSize, "", "PLANAR")
euc_eop = r"in_memory\eucdist_eop"
eopslopes = arcpy.Raster(masked_dem3) - (float(SinkVal) / (arcpy.Raster(euc_eop) + 1))
eopslopes.save(r"in_memory\eopslopes")

the above yeilds this error when run in teh script tool GUI from Caralog

PYTHON ERRORS:
Traceback info:
    File "script.py", line 428, in <module>
    eopslopes = arcpy.Raster(masked_dem3) - (float(SinkVal) / (arcpy.Raster(euc_eop) + 1))

Error Info:
ERROR 000824: The tool is not licensed.

Other attempts attempt (fails):

I've tried several things, that all fail. here's an example of one of those failed changes I tried to save the input rasters as raster layers :

    arcpy.MakeRasterLayer_management(masked_dem3, "masked_dem3_lyr", where_clause="", envelope="", band_index="")
    arcpy.MakeRasterLayer_management(euc_eop, "euc_eop_lyr", where_clause="", envelope="", band_index="")
    expression = r"arcpy.Raster(masked_dem3_lyr) - (float(SinkVal) / (arcpy.Raster(euc_eop_lyr + 1))"
    arcpy.gp.RasterCalculator_sa(expression, "eopslopes")

What is the correct formatting for arcpy.gp.RasterCalculator_sa(expression, "outras")?

.

.

.

In response to Jae's answer which I guess directly addressed my final bolded question above:

OK so arcpy.gp.RasterCalculator_sa() is out...but more to the spirit of my question, what then should one use in this case? my first attempt above works in IDLE, what formatting needs to be done to have it work in a script tool?


texts:
I want to process raster data (DEMs) using various algorithms of my own design.  I want to load (e.g.) a GeoTIFF, run an algorithm against its data, and save the results to a new GeoTIFF.  I would prefer working in C or C++ since they're fast and I'm familiar with them.  What would be the best way to go about doing this?

Here's some toy pseudo-code showing the kind of thing I want to do:

geotiff src = new geotiff("source.tif");
geotiff dst = src.clone();

for (int y = 0; y < src.height; y++) {
    for (int x = 0; x < src.width; x++) {
        if (
             (x < src.width - 1) &&
             ( src.data[x][y]   != src.NODATA ) &&
             ( src.data[x+1][y] != src.NODATA )
           ) {
            dst.data[x][y] = src.data[x][y] - src.data[x+1][y];
        } else {
            dst.data[x][y] = dst.NODATA;
        }
    }
}

dst.write("output.tif");


texts: Seeking Shapefiles of Restaurants in South Florida?

texts: I am doing a project on food trucks for my Intro to GIS class and am looking for locations of places to eat here in South Florida to overlay with places of business. I would like to make a map of which employees have no place to eat for lunch to provide to a non-profit organization that has a job-training program on a food truck.


texts: Not able to install postgresql-10-postgis-2.4 in Ubuntu 18.04

texts:
I am not able to install postgis.

$ sudo apt-get install postgresql-10-postgis-2.4 libgdal20 libhdf5-100
libnetcdf13 libgfortran4
Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies:
 libgfortran4 : Depends: gcc-7-base (= 7.3.0-16ubuntu3) but
 7.3.0-27ubuntu1~18.04 is to be installed
E: Unable to correct problems, you have held broken packages.

This is what it says finally.


texts:
Yesterday i downloaded an Aster L1B image from USGS's website. If i load it in QGIS, it doesn't get displayed properly(location wise). Even after i've enabled "On the fly" option, it's not reading the WGS84 coordinate system (Figured out this type of glitch happened in ArcGIS also  - How to import and project ASTER 1B data into Arcmap?). I tried to assign projection to the image but it didn't work (Got the message - hdf has no raster bands).

I tried the same in GRASS(r.in.gdal) and it got "Error: Selected band (1) doesn't exist"

But if i open the same image in Global Mapper, the image gets displayed perfectly.

How to solve this ?


texts: gdal_translate does not crop the GeoTIFF file

texts:
When I use gdal_translate to crop the GeoTiff file, the output file is saving the same as the input file. Nothing changed. The size needs to decrease at the boundaries I specify. What can I do?

from osgeo import gdal
import os

os.system('gdal_translate -projwin 3046477.80 4645536.98 3046546.08 4645571.24 input.tif output.tif')


texts:
I have created images in google earth engine , when I work in the code editor platform it works and it has the right polygon shape  but whenever I export it to my drive and download it, I don't get the full Image.

This is how it looks on GEE -

and this is what I get afterr exporting to drive -

this is the code I wrote;

var diff1=(period1.subtract(period2)).divide(period1);
var diff2=(period2.subtract(period3)).divide(period2);
var diff3=(period3.subtract(period4)).divide(period3);
var diff4=(period4.subtract(period5)).divide(period4);
var meandiff=(diff1.add(diff1).add(diff3).add(diff4)).divide(4);

var visDiff={min: 0, max: 1, palette: [
    'FFFFFF', '#ffcccc', '#ff8080', '#ff1a1a', '#e60000'
  ]};

Map.centerObject(period1,7);
Map.addLayer(diff1,visDiff,'Diff periods1-2');
Map.addLayer(diff2,visDiff,'Diff periods2-3');
Map.addLayer(diff3,visDiff,'Diff periods3-4');
Map.addLayer(diff4,visDiff,'Diff periods4-5');
Map.addLayer(meandiff,visDiff,'18Years');

 var toexport = diff3.visualize(visDiff)

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: toexport,
  description: 'diff3',
  scale: 500,

});

any idea what can be the problem?


texts: SQL queries, it is better to do it in QGIS or PostGIS

texts:
With QGIS I can manipulate data stored on a PostGIS Geodatabase.

In order to query to the geodatabase, is it better to do it from QGIS or directly in PostGIS?

What are the benefits of querying in PostGIS?


texts: Opening Raster with "bad" vertical coordinate system crashes program

texts:
I have an elevation TIFF file that comes from another program that seems to have something wrong with its Vertical Coordinate System. The coordinate system is EGM96 Geoid. I cannot open this TIFF in code, but I can open this TIFF in ArcMap. And in ArcMap I can see that the Vertical Datum is "EGM96 geoid".

If I select the predefined EGM96 Geoid coordinate system in ArcMap the Vertical Datum is "EGM96_Geoid". So the datum is the only difference I can see between the two within ArcMap. Also in ArcMap I can use Define Projection to update the TIFF's projection to the predefined EGM96 coordinate system and then I am able to open the TIFF in code without any issues. I do not want to have to run define projection on these TIFF's.

When I try and open the original TIFF in code the program crashes instantly without throwing any exceptions or logging any errors. Here is the code I use to open the TIFF.

IWorkspaceFactory raster_factory = new RasterWorkspaceFactory();
IRasterWorkspace raster_ws = raster_factory.OpenFromFile(ws_path, 0) as IRasterWorkspace;
IRasterDataset raster_ds = raster_ws.OpenRasterDataset("raster_name");
IRaster2 raster = (raster_ds as IRasterDataset3).CreateFullRaster() as IRaster2;

Is there some other way I should be handling this TIFF? Surely if the ArcMap application can open this Tiff, I should also be able to open it in C#.

I am using ArcMap 10.5 and Windows 10.


texts: How to speed up PostGIS clipping of out-db rasters?

texts:
I'm clipping (and unioning) thousands of tiled raster files that are stored in an out-db PostgreSQL/PostGIS database located on a SSD.  I want to use PostgreSQL but I found that PostGIS is taking almost three times as long to clip than the gdalwarp.exe utility.
Is there something I can do to get PostgreSQL to perform closer to gdalwarp.exe?
Using gdalwarp.exe to clip raster takes 0.36 seconds versus the SQL query that takes about 1.0 second to complete. The envelope intersects with the raster. The raster is a 1024x1024 JP2OpenJPEG file.
I'm using Postgres 12.4 and PostGIS 3.0 on Windows 10.
I've tried adjusting various configuration settings in the Postgres.conf file (shared_buffers, work_mem, etc.) but couldn't find anything that made a difference.
I tried smaller tile sizes when setting up the out-db images but that actually slowed down the query.
gdalwarp.exe -crop_to_cutline -of GTiff -te 103.878150, 1.819093, 103.904514, 1.846015  rasterInput.jp2  clippedRaster.tiff

SELECT ST_Clip(rast, ST_MakeEnvelope(103.878150, 1.819093, 103.904514, 1.846015, 4326), true)
FROM _imagery
WHERE rid = 47300       -- rasterInput.jp2

EXPLAIN ANALYZE:
Index Scan using _imagery_pkey on _imagery  (cost=0.42..2.89 rows=1 width=32) (actual time=1010.391..1010.887 rows=1 loops=1)
Index Cond: (rid = 47300)
Planning Time: 0.414 ms
Execution Time: 1010.939 ms


texts: Is there a way to automatically populate attribute values in Qgis?

texts:
We are using a Qgis/PostGis combo, where Qgis is used for adding features with attributes into a PostGis database.
Is there a way to automatically fill the attribute data values with some project-specific constant values (Qgis project name, user name)?


texts: Merge layer files in folders and subfolders with arcpy

texts:
I work with arcview 10.3. I try to merge layer files in order to get one layer file that include all the features in the other layer files. All layer files as different names, and they all spread in big folder that divided to a lot of sub folders. My code is:

import arcpy,os,sys,string,fnmatch
import arcpy.mapping
from arcpy import env

rootPath = r"C:\Project\layers"
pattern = '*.lyr'
lyr2merge = []
counter = 0
for root, dirs, files in os.walk(rootPath):
    for filename in fnmatch.filter(files, pattern):
        lyr2merge.append(os.path.join(root, filename))
        counter = counter + 1
 arcpy.Merge_management(lyr2merge, r"C:\Project\layers\layer_total.lyr")
print 'merge'
print counter

i got an error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Table View.
ERROR 000840: The value is not a Table View.
ERROR 000840: The value is not a Table View.
Failed to execute (Merge).

i can merge shapefiles Merge shapefiles in folders and subfolders with arcpy but after that i will have to save it as layer file.
Is it possible?


texts: How to enable "Label every part of multi-features" with pyQGIS?

texts:
Regarding the 'Label every part of multi-features' check box located in the 'Labels // Rendering // Feature Options' layer properties UI...

Can this be done with pyQGIS?

These two commands enable the labeling, and allow color/font/etc. to work fine.

self.contourLayer.setCustomProperty("labeling", "pal")
self.contourLayer.setCustomProperty("labeling/enabled", "true")

Tests with these three commands do not enable multi-feature labels.

self.contourLayer.setCustomProperty("labeling/displayAll", "true")
self.contourLayer.setCustomProperty("labeling/labelPerPart", "true")
self.contourLayer.setCustomProperty("labeling/limitNumLabels","false")


texts: How do I join an area that has been split into many .dwg or .gml files?

texts: How do I merge a bunch of adjacent (1) dwg or (2) gml files into a continuous map using free software or online? Thanks.


texts:
I have a TopoJSON file which describes electoral districts in Montreal, and I would like to produce an adjacency matrix which would tell me which districts share a boundary with which other districts. Is there an easy way to do this?

Edit: here is the file https://github.com/nicolaskruchten/mtlelection2013bydistrict/blob/master/districts.topojson


texts:
I am trying to use ArcMap's Field Calculator to populate values in a column based on a value in another column.  So, in pseudo-code:

if [field1] = 'Foo'
  set [field2] = 'Bar'
endif

Something like that.  Test for 'foo' in the first attribute, and set the record's second attribute to 'bar.'  Seems like it should be pretty easy to do, but I keep receiving a "user interrupt" error every time I try running anything.  I think I have a misunderstanding of how to build VBA expressions with the Field Calculator...


texts:
Is it possible to load a raster file (e.g. geotiff or erdas img) directly into a Leaflet map? If so, can you guide me to the correct documentation or method to do this? I am trying to do this with the imageOverlay() function in Leaflet currently.

var imageUrl = 'http://localhost/test/my_raster.tif',
        imageBounds = [[30.430224, -87.193133], [30.435867, -87.185944]];
        L.imageOverlay(imageUrl, imageBounds).addTo(map);

I am new to Leaflet.


texts: Creating Peirce quincuncial map in QGIS?

texts:
I would like to create a Peirce quincuncial projection of vector data using QGIS.

A similar question was asked here several years ago, without apparent satisfactory answer. Someone offered a link to Github code in a comment, but I don't have the coding knowledge to know how to implement it or whether/how it can be used with QGIS.

How can I use that link, or do this another way?


texts: Is there a version for Google Earth Pro x64 in Windows?

texts:
As the title says. I have looked more or less everywhere and couldn't find any definitive answer. If not for Windows, is there a x64 version for *nix (Ubuntu, Debian etc)?

Otherwise, with 32-bits Google Earth and 3 GB memory limit, is a SSD the best and only option to speed up things a little?


texts: Selecting multiple values with Select by Attributes in ArcGIS Desktop?

texts:
Can anyone help me with selecting multiple values from the attribute table?

I tried "Classes"='14'AND'07' but it doesn't work.


texts: Splitting polyline by point using ArcObjects?

texts:
I want to split gas pipe line by point and add two new feature and delete the old feature gas line.
I reach the code to split pipe line but does't work and stopped there
 ISet newFeaturesSetTemp = featureEdit.SplitWithUpdate(pPoint);

  public override void OnMouseDown(int Button, int Shift, int X, int Y)
{
    IMxDocument mxdoc = m_application.Document as IMxDocument;
    IActiveView activeView = mxdoc.ActiveView;
    IPoint pPoint = activeView.ScreenDisplay.DisplayTransformation.ToMapPoint(X, Y);
    IFeatureLayer pFlayer = null;
    IMap map = mxdoc.FocusMap;
    for (int i = 0; i < map.LayerCount - 2; i++)
    {
        ILayer layer = map.get_Layer(i);
        if (/*/layer is IFeatureLayer &&/*/ layer.Name == "Double_Endcap")
        {
            pFlayer = layer as IFeatureLayer;
            break;
        }
    }
    ///////////////////////////////////////////////////////////////////////////////////////////////////////
    IFeatureLayer flayer = null;
    for (int i = 0; i < map.LayerCount - 2; i++)
    {
        ILayer layer = map.get_Layer(i);
        if (/*layer is IFeatureLayer &&*/ layer.Name == "PePipe")
        {
            flayer = layer as IFeatureLayer;
            break;
        }
    }
    ITopologicalOperator topologicalOperator = pPoint as ITopologicalOperator;
    IGeometry geometry = topologicalOperator.Buffer(1);//(0.5);
    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.Geometry = geometry;
    spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
    //IFeatureSelection featureSelection = flayer as IFeatureSelection;//
    //featureSelection.SelectFeatures(spatialFilter, esriSelectionResultEnum.esriSelectionResultNew, true);//
    IFeatureClass featureClass = flayer.FeatureClass;
    IFeatureCursor featurecursor = featureClass.Search(spatialFilter, false);
    IFeature feature = featurecursor.NextFeature();
    ///////////////////////////////////////////////////////
    int iOID = feature.OID;
    ////Split the feature. Use the IFeatureEdit::Split method which ensures
    ////the attributes are correctly dealt with.
    IFeatureEdit2 featureEdit = feature as IFeatureEdit2;
    ////Set to hold the new features that are created by the Split.
    ISet newFeaturesSet = new SetClass();
    ISet newFeaturesSetTemp = featureEdit.SplitWithUpdate(pPoint);
    if (newFeaturesSetTemp != null)
    {
        newFeaturesSetTemp.Reset();
        for (int featureCount = 0; featureCount < newFeaturesSetTemp.Count; featureCount++)
        {
            IFeature newFeature = newFeaturesSetTemp.Next() as IFeature;
            if (newFeature.OID == iOID)
            {
                newFeaturesSet.Add(newFeature);
            }
        }
    }
    newFeaturesSet.Add(feature);
    ////////////////////////////////////////////////////
    if (feature == null)
    {
        MessageBox.Show("Select Pipe");
    }
    else
    {
        IFeatureClass pFeaturecalass = pFlayer.FeatureClass as IFeatureClass;
        IDataset pDataset = pFeaturecalass as IDataset;
        IGeoDataset pGeoDataset = pFeaturecalass IGeoDataset;
        ISpatialReference sr = pGeoDataset.SpatialReference;
        sr.SetDomain(-1000000, 40000000, -1000000, 40000000);
        IWorkspace pWorkspace = pDataset.Workspace as IWorkspace;
        IWorkspaceEdit pWorkspaceEdit = pWorkspace as IWorkspaceEdit;
        pWorkspaceEdit.StartEditing(true);
        IFeature pNewFeature = pFeaturecalass .CreateFeature() as IFeature;
        pNewFeature.Shape = pPoint;
        pNewFeature.Shape.SpatialReference = sr;
        pNewFeature.Store();
        pWorkspaceEdit.StopEditOperation();
        pWorkspaceEdit.StopEditing(true);
        activeView.Refresh();
    }
}


texts:
I'm trying to set up a loop that includes using items in the array along with string text to populate the where clause.  I'm getting the invalid sql expression for the clause.  Here's what I have so far.

# Flight Levels array

 flightLevelsArray = ["FL290", "FL300", "FL310", "FL320", "FL330", "FL340", "FL350", "FL360", "FL370", "FL380", "FL390", "FL400", "FL410", "FL420"]

queryFL290 = "FL IN ('FL285', 'FL286', 'FL287', 'FL288', 'FL289', 'FL290', 'FL291', 'FL292', 'FL293', 'FL294')"
queryFL300 = "FL IN ('FL295', 'FL296', 'FL297', 'FL298', 'FL299', 'FL300', 'FL301', 'FL302', 'FL303', 'FL304')"

for flightLevel in flightLevelsArray:
flightLevelField = "!"+flightLevel+"!"
rangeToQuery = "query"+flightLevel

arcpy.SelectLayerByAttribute_management(out_layer_file, "ADD_TO_SELECTION", rangeToQuery)


texts: How to draw line pointing in compass bearing direction from point shapefile?

texts:
I captured waypoint data in the field including the compass azimuth when looking at some summits.

I want to use these waypoints, that I have now in Shapefile format, to draw line shapefiles that point in the compass azimuth that I defined for each waypoint.

Do you know an elegant solution for this problem?


texts: GPServer does not understand %VARIABLE% sintax

texts:
I'm creating a model in order to publish it as a geoprocessing service. The model runs successfull in ArcMap but it crashes when I run the published GPService.

The model has a parameter called FUSO, which can be equals to 23 or 24. This parameter selects for instance which geodatabase table I'll need to use as a route source.

The GPServer errror is Error executing tool.: The value does not exist.\nERROR 000732: Input Features: Dataset GEO.VWG_LOGRADOURO_LIN_%FUSO% does not exist or is not supported\nWARNING 000970: The join field id_logradouro_%fuso% in the join table TABELA_DISTANCIAS is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.

So I guess that despite of ArcMap be able to interprets %FUSO% as the value of the parameter FUSO, the GPServer does not.

Refer to images bellow to see how I built the model.

How can I pass throught it?


texts: What can cause an ArcIMS server process to run away with the CPU?

texts:
We have an old application running in ArcIMS (don't ask, we're getting rid of it as soon as we can :)). All who had anything to do with the application's design or programming are long gone. We're trying to make this 9.3.1 application run in version 10 so we can upgrade our production servers to ArcGIS Server 10. Our ArcIMS runs on a server along with ArcGIS Server.

Every 24 hours or so one or more arcimsserver.exe processes go rogue and take up 99 - 100% of the CPU. On the ESRI forums this type of question doesn't seem to get a response, and ESRI support has failed to help is for the last several incidents. Does anybody here have a suggestion for what to check?

EDIT: The data includes image and map services. The data is in shapefiles and SDE.

EDIT: After some testing and exploring, we found that a database server had been moved since the application was built. The server connection was buried in the application. When we corrected it and restarted, the problem did not occur as often. After some more digging, we found other obsolete data connections. We have corrected them and now the problem occurs only about once every ten days or so.


texts: Algorithm to simplify a Polygon/Series of Points to a given number of Points

texts:
I have a polygon (or let's call it just a series of points). I want to simplify it (reduce the amount of points) to a given number of points.

Currently I use Douglas Peucker and try to guess the correct epsilon. That works more or less (I get approximately the given amount of points). The biggest problem is, that I have to call DP over and over again. Since I want to massively reduce the points, e.G. from 10.000s to 30 the process is a bit frustrating.

So I thought, that I might create my own version of DP, were I stop when I reach the amount of points I am looking for. Are there any existing or maybe better solutions available?


texts: Changing Field Name/Apply Symbology from layer using ArcPy?

texts:
Using ArcMap 10.1, I want to use the Apply Symbology From Layer tool, however my field names are not the same. I'm having trouble understanding how to change the Field Names in a layer.

I tried the Make Features Layer tool and it did not rename my variables. My variable names are basically the "LayerName"+"Variable Name".

Has anyone else had this problem?


texts:
I get the following error when trying to convert a .las file using las2txt with the LASTools plugin:
My system is set to UTF-8.

QGIS code revision: 4470baa1a3
Qt version: 5.11.2
GDAL version: 2.4.1
GEOS version: 3.7.2-CAPI-1.11.0 b55d2125
PROJ version: Rel. 5.2.0, September 15th, 2018
Processing algorithm…
Algorithm 'las2txt' starting…
Input parameters:
{ 'ADDITIONAL_OPTIONS' : '', 'CPU64' : False, 'GUI' : False, 'INPUT_LASLAZ' : '\\\\scansrv.scanergy.local\\Folder Redirection\\ia\\My Documents\\Svavelsö GIS\\laserdata\\10C031_65925_6825_25.las', 'OUTPUT_GENERIC' : 'TEMPORARY_OUTPUT', 'PARSE' : 'xyz', 'VERBOSE' : False }

LAStools command line
C:\LAStools\bin\las2txt -i "\\scansrv.scanergy.local\Folder Redirection\ia\My Documents\Svavelsö GIS\laserdata\10C031_65925_6825_25.las" -o "C:/Users/ia/AppData/Local/Temp/processing_10592c2961fb4f2b8e6471b587b96e4a/1f021684e02c4a05859479860806ffba/OUTPUT_GENERIC.file"
LAStools console output
Traceback (most recent call last):
File "C:/Users/ia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\LAStools\LAStools\las2txt.py", line 55, in processAlgorithm
LAStoolsUtils.runLAStools(commands, feedback)
File "C:/Users/ia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\LAStools\LAStoolsUtils.py", line 61, in runLAStools
feedback.pushConsoleInfo(output.decode("utf-8"))
UnicodeDecodeError: 'utf-8' codec can't decode byte 0x86 in position 199: invalid start byte

Execution failed after 0.15 seconds

Loading resulting layers
Algorithm 'las2txt' finished


texts: How to change the raster NoData value with ArcObjects?

texts:
I am trying to change the a raster's NoData value with ArcObjects and keep coming up short.

Here is the code that I thought would work:

IRasterDataset _dataset = ...get the dataset...
object[] noDataValue = ...set the no data values by raster band...
var rasterbandCollection = (IRasterBandCollection) _dataset;
for(int i = 0; i < rasterbandCollection.Count; i++)
{
    ((IRasterProps) rasterbandCollection.Item(i)).NoDataValue = noDataValue[i];
}

This value persists during the ArcMap session, but does not save the changes to disk it appears. Any ideas on how to make this permanent (without having to save a second raster dataset)?


texts:
In Geoserver, a have a SLD style with a rule that looks like this

    <Rule>
        <ogc:Filter>
            <ogc:PropertyIsGreaterThan>
                <ogc:Function name="env">
                  <ogc:PropertyName>
                      property
                  </ogc:PropertyName>
                </ogc:Function>

                <ogc:Function name="env">
                  <ogc:Literal>value1</ogc:Literal>
                  <ogc:Literal>0</ogc:Literal>
                </ogc:Function>

            </ogc:PropertyIsGreaterThan>
        </ogc:Filter>
        <PolygonSymbolizer>
            <Fill>
                <CssParameter name="fill">
                    #ff0000
                </CssParameter>
                <CssParameter name="fill-opacity">
                    0.5
                </CssParameter>
            </Fill>
        </PolygonSymbolizer>
    </Rule>

In the WMS request, i send an env parameter, specifing the property and the value1 parameters for dynamic sld substitution. The property parameter gets a value of an attribute name that exists on the layer I want to display the style with.

My solution won't work. I get an exception, saying:

  The requested style cannot be used with this layer. The style
  specifies an attribute property and the layer is: ...

Is there a way to set the attribute name dynamically via sld substitution in Geoserver? And if yes, how can I achive that?

I am using Geoserver 2.7.2

Btw. this works:

            <ogc:PropertyIsGreaterThan>

                  <ogc:PropertyName>
                      Einwohner
                  </ogc:PropertyName>

                <ogc:Function name="env">
                  <ogc:Literal>value1</ogc:Literal>
                  <ogc:Literal>0</ogc:Literal>
                </ogc:Function>

            </ogc:PropertyIsGreaterThan>


texts: Creating new minimum and maximum values to be consistent across 4 DEMs

texts:
I have 4 DEMs (LIDAR, 5M, 10M AND 30M resolutions) all with the same geographical extent and I want to compare the differences in resolution using 4 layouts. I want to have the same minimum and maximum elevations for each DEM. They vary slightly but all the minimum values are around 1520 and all the maximum values are around 3450. Someone told me I can use the Con tool to make all elevations below and above certain values NoData but I am having difficulty figuring out how to properly use this tool (my SQL expression is always giving me an error when I verify it so I am not even sure how to go about doing this). I suppose I could just change the numbers in the legend but I want to do it properly and actually have the same range of elevation values for each DEM.

How should I proceed?

I am using Esri's ArcMap 10.7.1.


texts: Land Distance from many coastal points to one point following the land route

texts:
I have a lot of coastal points around Mediterranean sea and One point in South Africa.I want to calculate the distance from all of these Mediterranean points to the one in South Africa but only following the land.However, the distance tools measure this distance passing from the sea(for example for the south France points).

How I can measure this land distance route?

I have been advised to use the path distance tool,but it seems quite difficult for me.

May you explain me the steps or could I have another recommendation?


texts:
is there anyway to set tolerance in Postgis, in Arcgis tolerance is the minimum distance which could separate two points, if the distance between them is smaller than the tolerance so the points are considered overlapped. i need to set this parameter because i get some weird results using functions like ST_Intersects, ST_equals, ST_difference. in this conditions i check intersecting points by ST_Distance(p1 , p2) < tolerance but i would like to get a more robust and general solution.

in my use case i need only 1cm to 1mm precision, but for example sometimes ST_Intersects do not give expected results because of some microns (0.000001).

SELECT ST_Intersects('POINT(0.00000000000000001 0)'::geometry, 'LINESTRING ( 0 0, 0 2 )'::geometry);


texts: Min NDVI estimation from NDVI time series

texts:
I'm trying to estimate max min NDVI values from the NDVI datasets that came from Landsat images 4-5-7-8. For the max NDVI I used the "qualityMosaic" and it seems that it works okay, but for the min NDVI value I get an error message

  MIN NDVI: Layer error: Computation timed out

I've attached the GEE code.

// Function to cloud mask Landsat 8.
var maskL8SR = function(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = ee.Number(2).pow(3).int();
  var cloudsBitMask = ee.Number(2).pow(5).int();
  // Get the QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(
            qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image
      // Scale the data to reflectance and temperature.
      .select(['B4', 'B5'], ['RED', 'NIR']).multiply(0.0001)
      .addBands(image.select(['B11'], ['Thermal']).multiply(0.1))
      .updateMask(mask);
};

// Function to cloud mask Landsats 5-7
var maskL57SR = function(image) {
  var qa = image.select('pixel_qa');
  // Second bit must be zero, meaning none to low cloud confidence.
  var mask1 = qa.bitwiseAnd(ee.Number(2).pow(7).int()).eq(0).and(
      qa.bitwiseAnd(ee.Number(2).pow(3).int()).lte(0)); // cloud shadow
  // This gets rid of irritating fixed-pattern noise at the edge of the images.
  var mask2 = image.select('B.*').gt(0).reduce('min');
  return image
      .select(['B3', 'B4'], ['RED', 'NIR']).multiply(0.0001)
      .addBands(image.select(['B6'], ['Thermal']).multiply(0.1))
      .updateMask(mask1.and(mask2));
};

// find all data and filter them by date
var lst5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
    .filterDate('1984-10-01', '2011-10-01')
    .filterBounds(roi)
    .map(maskL57SR)
    .map(function(image){return image.clip(roi)});

var lst7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
  .filterDate('2011-10-01', '2013-04-07')
  .filterBounds(roi)
  .map(maskL57SR)
  .map(function(image){return image.clip(roi)});

 var lst8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
     .filterDate('2013-04-07', '2018-05-01')
    .filterBounds(roi)
  .map(maskL8SR)
    .map(function(image){return image.clip(roi)});

// Combine all landsat data, 1985 through 2015
var L4578 = ee.ImageCollection(lst5.merge(lst8));
L4578 = L4578.merge(lst7);

// create function to add NDVI using NIR (B5) and the red band (B4)
var getNDVI = function(img1){
  return img1.addBands(img1.normalizedDifference(['NIR','RED']).rename('NDVI'));
};

// map over image collection
var l8ndvi = L4578.map(getNDVI).select('NDVI');

///////////////////////////////////////////////////////////////
// for each pixel, select the "best" set of bands from available images
// based on the maximum NDVI/greenness
var composite = l8ndvi.qualityMosaic('NDVI').clip(roi);

print (composite);

//set visualization parameters for Maximum NDVI

var ndviParams = {min: 0, max: 1, palette: ['red', 'yellow', 'green']};
Map.addLayer(composite.select('NDVI'), ndviParams, 'Maximun NDVI');

///////////////////////////////////////////////////////////////////////
//For each pixel, select the "worst" set of bands from available images
//based on minimum NDVI
// reduce the image collection to one image by taking MIN of the rasters
var reducemin_NDVI = l8ndvi.reduce(ee.Reducer.min());

// get min NDVI values by ROI polygon

var NDVImin = reducemin_NDVI.reduceRegions({
  collection: roi,
  reducer: ee.Reducer.min(),
  scale: 30 // the resolution of the dataset

});

//set visualization parameters for MIN NDVI

var ndviParams = {min: 0, max: 1, palette: ['red', 'yellow', 'green']};
Map.addLayer(NDVImin, ndviParams, 'MIN NDVI');


texts: How to perform analysis on point datasets in ArcMap?

texts:
I'm new with ArcGIS and currently I'm working on a project involving CAFO waste storage type, land cover, and Nitrate levels in Northern Missouri. I have land cover as a raster data set, with six different values indicating terrain. Nitrate level is represented by point data split into ranges. CAFO type is split into 3 different values. I'm trying to "perform analysis" on CAFO type with Nitrate level. They are both point data sets and I'm really just trying to relate them. Any suggestions?

Nitrate (mg/L):
0-.25, .30-1.0, 1.1-1.9, 2.0-9.9, >10

CAFOs: Storage Lagoon, Underfloor Pits, No Storage

Main attributes for land cover are forest, pasture, and open water.

I'm trying to see if there's any correlation between higher Nitrate levels and a specific storage type (for example, if Storage Lagoons are associated with higher Nitrate levels). Is there any analysis function that would provide me with this information? I could also relate land cover to Nitrate levels, but I think it might be more complicated stepping outside of the point-to-point.

Here's a photo of what I've got so far, I know, it's a mess...


texts: Delete isolated polygons smaller than a specified area

texts:
I want to delete isolated polygons (only polygons not adjacent to other polygons) smaller than a specified area. Like the polygons displayed here in the white space.

I can do a select analysis - area less than, but that removes ALL small polygons. The small isolated polygons are not that important as they represent tiny values.

I have an ArcEditor license.


texts: Read R-Tree from PostGIS enabled database

texts: There is a method available to get the largest rectangle from the R-Tree (see this question). Are there methods available that allow one to read the entire tree?


texts: How to sum the area depending from another field value in field calculator?

texts:
I have a landcover-layer with an attribute of the administrative region (GEN). I am looking for an expression in the field calculator to calculate the summarized area depending on the GEN-field value.
I know I could do it with dissolve tool or with group stats, but I would like to have the GEN area sum as a field value in my attribute table.


texts: Data Insertion from .CSV file to a Feature Class attribute Table in ArcGIS for Desktop?

texts:
I am working on water network in ArcGIS 10.2. I am trying to insert data from a data.csv file to an existing feature class House Meter in personal GDB.

One method I know to use is the join data.csv file on the basis of common field with House meters.

After that use field calculator for data insertion. This method takes too much time because I have to insert data in almost 15 fields from CSV file.

I want to insert data without join data.csv file and quickly.

Does anyone know of a script or expression for this query?


texts:
I'm in arcpy scripting for couple of years and after this time i've got a whole lot of scripts, packages and so on. The problems is that many of this scripts should be used over and over, many of them use the same logic (couple of point-to-point spatial joins with different limitations). Additional problem is that I reuse some part of this scripts: extentded arcpy.da.Cursor class for example (not only this one, it's about 5 classes and 5 modules with functions. So, 10 "scripts" in total)

Now the time to orginize this mess has came. At this moment I use project structure. One project - one toolset with it's own pack of scripts. I just copy from base location to project's one if I need my class in another project

Additional problem is that this "package" should be portable as easy as possible.

Almost all of this scripts/tools are used inside ArcMap environment.

Now I think of a structure like this:

_docs                              # sphinx-generated docs
gptools                            # main "toolbox"
    --- __init__.py                # access everything from outside and generate docs
    --- gptooset1                  # same logic script. Just as ArcMap toolboxes.
        --- __init__.py            # i could import tools from here
        --- gptooset1.tool1.py     # single tool
        --- gptooset1.tool2.py
    --- gptooset2
    --- ...
scripts                            # main wideused classes, functions, etc.
    --- __init__.py                # use it from outside
    --- class1.py
    --- ...
    --- functions1.py
    --- ...
site-packages                      # third-party modules, should be imported from here
    --- module1
toolboxes                          # python-toolboxes to access tools from ArcMap

There is few issues. First of all I can't import class1 to gptoolset1.tool1.py with relative import. There is seems to be only one way to do it through sys.path.append() which is

The second, the main one is that this structure is against Python basic structure.

scripts
    --- class1.py
    --- functions1.py
tool1.py

# tool1.py
from scripts import class1
from scripts import functions1

Does anyone has an experience of organizing such large packages? How to organize it well from Python point of view and keep it as simple is possible? (bulk of sys.path.appends would solve it but its awkward.)

UPD
I don't think it's a general software question because it's ArcGIS specific. Python "software" uses up->down (complicated->easy) logic, when main modules are uses "small pieces of logic" by importing it. It gives a strict scheme (packages, subpackages and so on).

ArcGIS uses reversed logic: tools (python modules actually) are on the bottom of project structure and oftenly uses one common parts of code. If you put tools on top of your project you would get a mess of scripts or should copy-paste common parts from one location to anoter.

The question is how to combine this two completely different approaches. Why not SO? Because it's not that frequent case in python scripting but it is a part of ArcGIS environment.


texts: Applying symbology from style file using ArcPy?

texts:
Every guide/tutorial I found explains how to apply symbology to a layer using the "arcpy.ApplySymbologyFromLayer_management". Instead, I do have a style file, and I want to match the fields on the attribute file with the icons on the style file, through python (point layer, unique values symbology).

Is it even possible to do that through python?

I ask that because there are over 100 icons on the style file, each one for a possible value. If I use Apply Symbology from Layer, I will end up with dozen empty icons on the legend, which I have to remove manually.  So if anyone also knows how to remove icons that the count is zero on the symbology (through python), it also works


texts: Importing OpenStreetMap historical data with ogr2ogr?

texts:
I am trying to import historical OpenStreetMap data obtained here:
http://osm.personalwerk.de/full-history-extracts/2014-11-24/continents/
I downloaded the Australia-Oceania file, I am trying to create a spatialite database for Papua New Guinea, the same way I normally do with current data:
ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat 139.73 0.37 157.64 -12.7 png-full-historic.db australia-oceania.osh.pbf

Instead of the expected outcome, I get the following error:

Error: unsupported required feature : HistoricalInformation
ERROR 1: An error occured during the parsing of data around byte 158

The same command works perfectly with current OSM data.
Can anyone help me to identify what I am doing wrong?
What is the mentioned feature "HistoricalInformation"?


texts: Conditional Reclassification using Raster Calculator

texts:
There are several questions about conditional reclassification on the Stack Exchange, but I haven't been able to find exactly what I need...

In my workspace, I have a list of raster layers (8 bit unsigned integer) that have values from 0-255. I am attempting to use conditional statements within the raster calculator to create a binary (0 and 1) layer in which only values between 90-100 are reclassified to 1 and the remaining values are reclassified to 0.

The code snippet I am working with can be seen below

# IMPORT MODULES
import arcpy, os

# SET GEO_PROCESSING ENVIRONMENTS
arcpy.env.workspace = r"C:\Users\BP_Mosaics"
outws = r"C:\Users\90BP_Mosaics"

# CREATE A LIST OF THE MOSAIC RASTER LAYERS
rasterList = arcpy.ListRasters("*", "tif")

# CREATE LOOP TO RECLASSIFY ALL BP MOSAICS INTO 0s and 1s
# CREATE LOOP TO RECLASSIFY ALL BP MOSAICS INTO 0s and 1s
for i in rasterList:
    inputRaster = (os.path.join(env.workspace, i))
    filename = i[-18:-4] + '_01.tif'
    outname = os.path.join(outws, filename)
    arcpy.gp.RasterCalculator_sa("Con(((" + inputRaster + ">= 90) & (" + inputRaster + "<= 100)), 1, 0)", outname)

I am not seeing any immediate errors in the code, however something does not seem to be working properly as I receive the following error message...

arcgisscripting.ExecuteError: ERROR 000539: (null)
Failed to execute (RasterCalculator).

Would anyone happen to know what is happening here?


texts:
I'm printing a series of maps using the Atlas feature in QGIS.  Everything is working fine, except when I print the atlas all of the elements are offset from where they are in the composer.
I can work around by printing the atlas, then gaming the position so the elements  show up where I want, but it takes a bit of monkeying that shouldn't be necessary.
If I use the regular "Export as PDF" button, I get a copy of the current atlas page with everything aligned appropriately.  Using the "export Atlas as images also produces correctly aligned images, but then need to be converted and assembled into a book.

But when I print, which I need to do to get the entire atlas, and not just the current page, all of the map elements appear shifted down and to the right by about 1/3 inch.

It's a little difficult to see the offset in the images because of the white page on white background, but if you look closely, the map elements are all shifted down and to the right, and as a consequence, the margins are not even.

By measuring the shift and adjusting the positions, I can get things close, but this is really clunky, and doesn't produce as nice of results as I can get by specifying positions and sizes of the elements.

How can I force the atlas to print like it's displayed?


texts:
I have several files (shapefiles). Table 1 has only one row with several fields whereas Table 2 has several rows.
What I want is to take a value from Table 1 and populate all rows in Table 2. I used a SearchCursor to find the value I want and an UpdateCursor to populate the rows.

import arcpy

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Ecodis_clip_split"
#set workspace to folder with values

rows = arcpy.SearchCursor ("Table1.shp", "", "", "area_ha", "")
#search for the value in attributetable

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Split_wetlandt"
#set workspace to folder with files which need to be updated

rows = UpdateCursor ("Table2.shp", "", "", "area_eco", "")
#Update all rows of file with value from SearchCursor above

Running this, I keep getting this Errormessage:

Traceback (most recent call last):
File "<module1>", line 17, in <module>
NameError: name 'UpdateCursor' is not defined

If something is not defined, this means that Python sees it as a variable?
Even if I could get rid of the errormessage: would it populate all the rows in Table 2 with the value or do I need to do something else?
Note that I just started to learn Python and never scripted anything what so ever...

Edit 6/17/14 11:32 am

@ Aaron, would this replace the values in Table 2 with my searched values from table 1?

import arcpy
arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Ecodis_clip_split"
#set workspace to folder with values

rows = arcpy.SearchCursor (table 1", "", "", "area_ha", "")
#search for the value in attributetable

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Split_wetlandt"
#set workspace to folder with files which need to be updated

with arcpy.UpdateCursor (Table 2 ["area_eco"])
#Update all rows of file with value from SearchCursor above
for row in cursor:
    area_eco.replace(SearchCursor)

EDIT 6/27614 12:04pm

I updated my script to this.
Somehow it deleted my whole attributetable in Table 2 but I got it restored. I run it again and now nothing happens. It does not update the rows, the value in the rows stays 0.

import arcpy
arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Ecodis_clip_split"
#set workspace to folder with values

rows = arcpy.SearchCursor ("Table1", "", "", "area_ha", "")
#search for the value in attributetable

arcpy.env.workspace = "C:\\Users\\Documents\\Scoring\\Wetland"
#set workspace to folder with files which need to be updated

rows = arcpy.UpdateCursor ("Table2", "", "", "area_ha", "")
del rows


texts: Getting number of objects in attribute table using PyQGIS?

texts:
I have a layer and I would like to access the number of rows the attribute table has from my code. Is there a way to do that without doing it like this:

Index = 0

for Feat in Layer.getFeatures():
    Index += 1

print Index


texts: Wrong Display in Line Symbol and Gradient Color

texts:
I create a line symbol with specific gradient color and below control template:

<ControlTemplate>
        <Path  x:Name='Element'   StrokeThickness='10' UseLayoutRounding='False' StrokeStartLineCap='Round' StrokeEndLineCap='Round' StrokeDashCap='Round' StrokeLineJoin='Round'  RenderTransformOrigin='0.5,0.5' StrokeDashArray='0' Stretch='Fill' >
            <Path.Stroke>
                <RadialGradientBrush GradientOrigin='8,0' Center='0,0' RadiusX='-4.5'>
                    <RadialGradientBrush.RelativeTransform>
                        <CompositeTransform CenterY='0.5' CenterX='0.5' />
                    </RadialGradientBrush.RelativeTransform>
                    <GradientStop Color='Transparent' Offset='0.588'/>
                    <GradientStop Color='Red' Offset='1'/>
                </RadialGradientBrush>
            </Path.Stroke>
        </Path>
</ControlTemplate>

Line in blend

and result is:

Line Symbol in map

What does reason of bad diplay in gradient color?


texts: Extract by Mask giving ERROR 999999 in ArcGIS Desktop?

texts:
In ArcGIS Desktop 9.3 I have 3 features I'm trying to use to extract by mask a DEM raster layer. After I extract by mask the layers visually don't line up. When I try to do further analyses (e.g., extract to point, slope, summary stats, etc.) it is clear that my mask does not line up with my features because I get -999999 as values in further analyses.

What went wrong? and what can I do to fix it?

I'll also note that my features were originally drawn polygons that I converted to features -- could this have anything to do with the problem?

STEPS

Create New layer
Draw Polygon
Manually adjust polygon length & width
Rotate polygon
Copy polygon 2x so that I have 3 identical polygons
Place polygons where they belong on map (using road features layer as guide)
Convert Graphics to features (using DEM as source data for coordinate system)

Selected all 3 drawn polygons at once and then selected to convert graphics to features

Add new features as layer

Agreed to the menu that popped up asking to add my converted features as a layer

Extract by Mask using new layer features as input and DEM for raster input

Properties of both Layers (I didn't manually define or adjust)

DEM Raster:

Spatial Reference: NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
linear Unit: Foot_US (0.304801)
Angular Unit: Degree (0.017453292519943295)
False_Easting: 2000000.002616666
False_Northing: 0
Central_Meridian: -79
Standard_Parallel_l: 34.33333333333334
Standard_Parallel_2: 36.16666666666666
Latitude_ Of_ Origin: 33.75
Datum: D_North_American_1983

Polygon Features:

Projected Coordinate System:     NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
Projection: Lambert_Conformal_Conic
False_Easting:  2000000.00261667
False_Northing: 0.00000000
Central_Meridian:   -79.00000000
Standard_Parallel_1:    34.33333333
Standard_Parallel_2:    36.16666667
Latitude_Of_Origin: 33.75000000
Linear Unit:    Foot_US

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Overall Data Frame properties:

NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
Projection: Lambert_Conformal_Conic
False_Easting: 2000000.002617
False_Northing: 0.000000
Central_Meridian: -79.000000
Standard_Parallel_1: 34.333333
Standard_Parallel_2: 36.166667
Latitude_Of_Origin: 33.750000
Linear Unit: Foot_US

GCS_North_American_1983
Datum: D_North_American_1983

Picture demonstrating alignment issue (2 of 3 drawn features shown)

Black/grey = DEM mask layer; Red outline = polygon feature outline; Blue = rivers; Green = roads


texts:
I am new DIVA - GIS user. I want to import country .grd file from diva - gis website. I tried with Bolivia, and it works fine, but it doesn't work with Bosnia (which actually is my subject of analysis), although the files are there.

Does anybody know what could be the issue and how to resolve it?


texts: Using MBTiles, GeoServer, and OpenMapTiles

texts:
I can’t understand what’s the matter:

I downloaded the map from https://openmaptiles.com/ in the mbtiles format

Installed GeoServer and extensions

Published layers

Several layers are displayed, the rest are displayed in white.

Why can't I get the rest of the layers?

Blockquote2019-12-31 11:09:44,274 INFO [geoserver.wms] -
Request: getMap
    Angle = 0.0
    BaseUrl = http://localhost:8080/geoserver/
    Bbox = SRSEnvelope[-2.998352410887321E7 : 2.998352410887321E7, -3.002261605686392E7 : 3.002261605686392E7]
    BgColor = java.awt.Color[r=255,g=255,b=255]
    Buffer = 0
    CQLFilter = null
    Crs = PROJCS["WGS 84 / Pseudo-Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","3857"]]
    Elevation = []
    Env = {}
    Exceptions = application/vnd.ogc.se_inimage
    FeatureId = null
    FeatureVersion = null
    Filter = null
    Format = image/png
    FormatOptions = {}
    Get = true
    Height = 768
    Interpolations = []
    Layers = [org.geoserver.wms.MapLayerInfo@4840559c]
    MaxFeatures = null
    Palette = null
    RawKvp = {EXCEPTIONS=application/vnd.ogc.se_inimage, FORMAT=image/png, REQUEST=GetMap, SRS=EPSG:3857, BBOX=-29983524.10887321,-30022616.05686392,29983524.10887321,30022616.05686392, VERSION=1.1.1, STYLES=, WIDTH=767, SERVICE=WMS, HEIGHT=768, TRANSPARENT=true, LAYERS=main_map:housenumber}
    RemoteOwsType = null
    RemoteOwsURL = null
    Request = GetMap
    RequestCharset = UTF-8
    ScaleMethod = null
    Sld = null
    SldBody = null
    SldVersion = null
    SortBy = null
    SortByArrays = null
    SRS = EPSG:3857
    StartIndex = null
    StyleBody = null
    StyleFormat = sld
    Styles = [StyleImpl[ name=generic]]
    StyleUrl = null
    StyleVersion = null
    Tiled = false
    TilesOrigin = null
    Time = []
    Transparent = true
    ValidateSchema = false
    Version = 1.1.1
    ViewParams = null
    Width = 767


texts: How to draw line by scratching on map using OpenLayers?

texts: Is it possible to draw line by scratching on the map using OpenLayers? I tried some examples, all the examples are drawing by clicking. I want to draw by scratching, just like drawing in normal life on the paper. Is it possible to implement by OpenLayers?


texts: QGIS Heatmap Plug: Filling a radius with one color

texts: I'm using the heat map plug in and I want to fill in my circles with just one color that corresponds with a specific color along a gradient. My tinkering has not been fruitful and every tutorial I've found has used gradient color. Suggestions?


texts: Automatically select the default feature in TemplatePicker ARCGIS Javascript

texts:
I am using this example to implement the editing on my map.

Editing Without Widget ARCGIS Javascript

The code employs a templatepicker which uses the layers on the map to show the features which can be selected and hence editing can be done on those features.

            templatePicker = new TemplatePicker({
                featureLayers: layers,
                rows: "auto",
                columns: 2,
                grouping: true,
                style: "height: auto; overflow: auto;"
            }, "templatePickerDiv");

However, you have to select a feature yourself before starting to edit it. After selection, the selected feature is fetched as

selectedTemplate = templatePicker.getSelected();

This selectedTemplate is then selected to put the points on the map.

selectedTemplate.featureLayer.applyEdits([newGraphic], null, null);

The problem is, I want to select the default or first feature in the templatePicker automatically as it loads. I tried finding if TemplatePicker has any method to select any feature but I couldnt. Is there any way that I could do that?


texts: Calculating length of polygon in ArcGIS Desktop?

texts:
I am using ArcGIS 10.2.1 Desktop with an Advanced license.

My company has electric lines that go through people's yards. The areas around these lines require trimming occasionally, which necessitates acquiring easements on some properties. See attached image for an example of a work area.

I am trying to automate a process to measure a polygon or polyline (the main electric line) that crosses multiple polygons (parcels) so we can request a specific length of easement when necessary. Our current process involves creating a new polyline layer and making a new segment for the frontage of each property so the length is automatically calculated. However, this is a very slow and tedious process, as we have tens of thousands of parcels with main line.

I've tried converting the parcels polygon to lines (Feature to Line tool) and labeling based on length. However, the parcels layer is from the county, and it's pretty messy, especially in rural/non-regular areas (not visible on my attached screenshot). There are a lot of extra vertices and line segments that get pulled into the new line layer which confound the length measurements with this method. I've tried simplifying polygons to reduce/remove these extras, but simplifying too much will alter the shape of the parcel and many extraneous portions still seem to transfer over.

I've also tried buffering around the main electric line to create a polygon, then using the Identity tool with the parcels layer to split it based on property divisions. My thought was to create polygons that spanned the front of peoples' properties, then measure the length of those new polygons (we aren't concerned about width here, since that varies based on local ordinances, so it's calculated later). This produced a length field, but its values aren't even close to the actual length of each polygon.

However, I've been unable to find a way to automate the measuring of polygon length. Minimum Bounding Geometry hasn't been helpful.

Is there a script I could use or a tool that I'm missing that could simplify this process of finding the length of a polygon?


texts: Is it possible to check the z factor of a hillshade in the layer properties in QGIS?

texts: I'm working with QGIS 3.18 and I want to compare the visual performance of different hillshade layers. For one that I already have I wanted to chek the building parameters, specially the z-factor but also the azimut, light height, etc. I thought I would find it under the info tab in the layer properties, but it looks that either I don't know where to look or that info is not available.


texts:
I have the following code but I can't seem to make it work, it only works when I only have a simple label.

from qgis.utils import iface

label_settings = QgsPalLayerSettings()
label_settings.drawBackground = True
fieldname_1 = "NAME_1"
fieldname_2 = "Data_cases"
label_settings.fieldName = '''concat('BASE: ' + to_string("''' + fieldname_1 + '''"),'\nL = ' + to_string("''' + fieldname_2 + '''") + 'm')'''
#label_settings.fieldName = ('''"NAME_1" || '\n '|| "Data_cases"''')
#label_settings.fieldName = "NAME_1"
label_settings.enabled = True
label_settings.placement= QgsPalLayerSettings.Center
label_settings.bufferDraw= True
label_settings.bufferSize= 1
#label_settings.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')

text_format = QgsTextFormat()

text_format.setFont(QFont("Arial", 12))
text_format.setSize(55)

background_color = QgsTextBackgroundSettings()
text_format.setColor(QColor("white"))

background_color.setEnabled(False)
text_format.setBackground(background_color )
label_settings.setFormat(text_format)

cases_layer.setLabeling(QgsVectorLayerSimpleLabeling(label_settings))
cases_layer.setLabelsEnabled(True)
cases_layer.triggerRepaint()


texts: Take mean of points within a radius using ArcGIS?

texts:
Is it possible to derive for every point in a shape file the spatial mean of all points within a radius of x km from that point in ArcGIS?

At the moment I can only find the point to raster tool which converts all points within a given pixel size to give the spatial mean, but I'd like to use circles instead of pixels.


texts: Cannot set a NULL value on the not null column applcnt_type_status with org.geotools.jdbc.JDBCDataStore

texts:
Using geotools to store shape files in a PostgreSQL database. Encountering an exception on the insert. Any suggestions?
 DataStore shapeDataStore = DataStoreFinder.getDataStore(file);
 String typeName = shapeDataStore.getTypeNames()[0];
 SimpleFeatureType inputType = shapeDataStore.getSchema(typeName);
 FeatureSource<SimpleFeatureType, SimpleFeature>
    source = shapeDataStore.getFeatureSource(typeName);

 FeatureCollection<SimpleFeatureType, SimpleFeature>
    inputFeatureCollection = source.getFeatures();

 Map<Object, Serializable> postGISparams = new HashMap();

 postGISparams.put("dbtype", "postgis");
 postGISparams.put("host", "localhost");
 postGISparams.put("port", 5437);
 postGISparams.put("database", "db");
 postGISparams.put("schema", "new_schema");
 postGISparams.put("user", "user");
 postGISparams.put("passwd", "password");

 // storage in PostGIS
 DataStore dataStore = DataStoreFinder.getDataStore(postGISparams);
   dataStore.createSchema(inputType);

 SimpleFeatureStore featureStore = (SimpleFeatureStore)
   dataStore.getFeatureSource(typeName);

 // write results
 featureStore.addFeatures(source.getFeatures());

 String typeName1 = dataStore.getTypeNames()[0];

 SimpleFeatureStore featureStore1 = (SimpleFeatureStore)
    dataStore.getFeatureSource(typeName1);

 featureStore1.addFeatures(inputFeatureCollection);

The exception:
java.io.IOException: Cannot set a NULL value on the not null column applcnt_type_status
 at org.geotools.jdbc.JDBCDataStore.insertSQL(JDBCDataStore.java:4165)
 at org.geotools.jdbc.JDBCDataStore.insertNonPS(JDBCDataStore.java:1903)
 at org.geotools.jdbc.JDBCDataStore.insert(JDBCDataStore.java:1731)
 at org.geotools.jdbc.JDBCInsertFeatureWriter.flush(JDBCInsertFeatureWriter.java:126)
 at org.geotools.jdbc.JDBCInsertFeatureWriter.write(JDBCInsertFeatureWriter.java:98)
 at org.geotools.data.InProcessLockingManager$1.write(InProcessLockingManager.java:296)
 at org.geotools.data.store.ContentFeatureStore.addFeature(ContentFeatureStore.java:302)
 at org.geotools.data.store.ContentFeatureStore.addFeatures(ContentFeatureStore.java:254)


texts: Geocoding with MMQGIS. How to use Google API Key?

texts:
I'm trying to geocode two csv files with 10 addresses each. When I tried with the first one, every address was located by Google, but when geocoding the second file I get this error:

Exceeded Daily Google Limit: 0 addresses geocoded

I got a Google API Key for geocoding and even after using this key, I still get the same error. These are the keys I got from Google

are these API keys the ones that I need?
If so, why does MMQGIS say that I've exceeded my daily limit?
Is there another way to geocode in QGIS?


texts: create a 3D line from 3D points in Arcgis

texts:
I have a large amount of 3D points that represent several (separate) linear elements. This means that these points partly have the similar/(same) x,y - coordinates and differ only in the height (z) (like power lines). But the single points have no relation to each other, i.e. the left point knows nothing about its right point or other points. Using the Near function I could only find the nearest point, but not the point on the other side of the point. The function Generate Near Table generates the information about both points, but unfortunately only works in 2D, so that superimposed line-like elements will mix in the 2D plane.
My goal now is to create a single line for each line-like element. In the beginning I tried to do this with the point2Line function, but the line was quite crisscross, so I tried to correct this with a small buffer. The next step would be to generate a centerline of this polygon, but ArcGIS does not seem to support the generation of a centerline from closed polygons.

Is there a way to create a polyline from 3D points with ArcGIS? A 2D line would be one approach, because then I could work it with this tutorial from next.


texts:
I want to transfer a SWMM Model to QGIS, the Subcatchments included.
I have tried:

Exporting the map as a .dxf File from SWMM and importing the file again in QGIS.
This only transfers the Nodes and the Links, not the Subcatchment Areas.
Installing the ImportEpaNetImpFiles extension in QGIS and importing the .inp file directly.
This leads to to Python Errors:

2021-03-30T14:11:41     WARNING    Traceback (most recent call last):
  File "C:/Users/Johann/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ImportEpanetInpFiles\main.py", line 131, in run
  epa2gis(filePath[0], self.epsg_crs)
  File "C:/Users/Johann/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ImportEpanetInpFiles\Epa2GIS.py", line 31, in epa2gis
  d.BinUpdateClass()
  File "C:/Users/Johann/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ImportEpanetInpFiles\readEpanetFile.py", line 13, in BinUpdateClass
  mm = getBinInfo()
  File "C:/Users/Johann/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ImportEpanetInpFiles\readEpanetFile.py", line 716, in getBinInfo
  nodeJunctionElevations.append(float(mm[1]))
 ValueError: could not convert string to float: '00180M01'

I am using QGIS 3.18 and EPA SWMM 5.1
I want to transfer the Subcatchemnt Data to be shown as areas in QGIS, so I could e.g. project the percentage of impervious area in a different color for each subcatchment. I am assuming there is a way to get the subcatchment data from swmm to polygon data in QGIS, to do this.
Is this possible?


texts:
There are a few poi that I can see on Josm that are not rendered to the final map. OSM database contains them at least for some month. I can' figure out why the doesn't appear on the map.

Example poi ids:

2022172177
2011325312
1579653536
1572586097
1750804585

Any idea?

Thanks,

Hubi


texts:
I have a couple of images that I have made with the service "maps made easy" but I have a problem. The images are projected with the CRS: "EPSG:4326, WSG 84"

My normal projection is "EPSG: 25833" so the images are slightly off in placement.

How can I reproject the images so they will fit?

I have tried to use this method, but I get an unreadable file back, I don't know what I am doing wrong.
How to reproject a raster in QGIS?


texts: How to stylize KML layers for Google Maps?

texts:
I am using QGIS 2.0 to convert my shapefiles into KML layers.

The purpose of converting my shapefiles into KML layers is to bypass the restriction imposed by Google Fusion Maps, which only allows one fusion table to be stylized.

However, I’m having issues retaining the layer style when exporting my shapefiles to KML layers. I’ve tried modifying the style tag found in the XML file, but nothing changes.

Does anyone have any suggestion on how to stylize KML layers?


texts: Summarizing raster values by vector category in QGIS (equivalent to ArcGIS Zonal Statistics as Table)

texts:
Is there a function in QGIS that will tabulate zonal statistics by vector category?
In ArcGIS, this can be done using Zonal Statistics as Table, but it looks like the QGIS Zonal Statistics function simply appends to the attribute table.  The question here QGIS equivalent to the ArcMap Zonal Statistics as Table is similar, but is asking about cell counts for raster values, whereas I am interested in average and standard deviation.


texts: CARTO Airship category widget returning NaN values

texts:
I am using Airship and carto-vl while working through this tutorial and the values in the widget are coming back as NaN. Here is my code:

<!DOCTYPE html>

  <link rel="stylesheet" href="https://libs.cartocdn.com/airship-style/v1.0.3/airship.css">

  <link rel="stylesheet" href="https://libs.cartocdn.com/airship-icons/v1.0.3/icons.css">

  <link href="https://libs.cartocdn.com/mapbox-gl/v0.48.0-carto1/mapbox-gl.css" rel="stylesheet" />

  <title>Airship Testing</title>

        <as-category-widget id="cities-widget" heading="Populated places" description="Maximun population of the most pouplated cities">
        </as-category-widget>

Here is what the app looks like:

Why are the values of the table not returning?

I thought there might be a bug in the tutorial, but I transferred the same methodologies over to a personal project and I got the same results.


texts:
Whats going on: I am working with postGIS and spatial datas. I am having trouble to fully understand the R-tree algorithm.

Example query: In which town is point P located?

My understanding so far is:
Each polygon has a boundingBox, this boundingBox has an idex. These boundingBoxes are then merged into another larger boundingBox which also has an idex. Continuing the process creates the R-tree. With the help of the R-tree the query is much faster. But ...

My question is at what point are these indexes created?
Does the R-tree exists per se?
Or is the R-tree created on the run time when I start a query?


texts: Mapping RBG color to coordinate possible?

texts:
I work with geologists who are trying to figure out how to use Qgis to batch apply colour values to coordinates so that we can visually see  differences or patterns.

From the attached photo of soil samples  I've created an Excel file containing all the sample numbers and their corresponding average RGB colour (using the Photoshop color picker). The geologist has all the samples and coordinates in Qgis but only knows how to manually apply an colour value to each sample point which will take foreverrrr.

Is there a way to do this automatically so that we can see a grid of slightly varying soil colours?

(Just as an example, and I hope this doesn't confuse the question, but it occurred to me that a pixel level view of a(ny) raster image represents essentially the same kind of effect we are looking for... except each pixel would be a real life soil sample colour, and the location would be the coordinates in Qgis where that sample/colour came from.  It doesn't have to be colour-up-against-colour like my example, but hopefully that helps illustrate the idea)

edit: I have tried doing this with an expression as advised by user PolyGeo but am getting a NULL error in the expression dialog after following these steps:  First we go to Layer Properties > Style Tab > Category ... then open the Expression Dialog > select Color dropdown > select color_rgb () option then input the variables into the brackets using the picker on the corresponding column labels Red, Green, and Blue. The resulting expression is color_rgb("Red", "Green", "Blue"). At the bottom of the dialog it says 'Exit Preview Error: NULL' This seems simple enough but we are clearly not doing something correctly.


texts: In Arcobjects, what is the difference between feature.shape and feature.ShapeCopy?

texts: What is the difference between Feature.shape and Feature.ShapeCopy in ArcObjects. I am using Feature.shape could you please tell me difference?


texts: GeoServer WMS / WMTS with Mapnik rendering backend

texts:
I am going to migrate our current TileCache WMS-T server to a current and up to date GeoServer / GeoWebCache product. We use Mapnik as our rendering back end and I would love to keep it this way. I run several layers written in Mapniks proprietary XML-Styling language... Therefor it would be wise to keep that.

Has anyone done this so far?
My idea so far is to run a lightweight WMS Server in front of Mapnik and then use that server to feed my GeoServers Tile storage....

Any other ideas?


texts: Invalid Output Workspace error when executing Slope_3d as WPS

texts:
I wrote a small toolbox that calculates the slope in degrees of a given .tif file, stores the ouput to a local folder on my computer and generates 2 by 2 tiles. The toolbox works fine, the publishing as WPS to ArcGis Server worked fine too. But when I try to execute the WPS in ArcMap several errors occur.

The first:

I intended to create a folder structure for the output tiles that every "zoomlevel" has it's own folder for better accessibility (folder 1 = 4 tiles, folder 2 = 16 tiles etc.):

import os, sys
import arcpy

####
## Environment Settings ArcPy
####

file_in = arcpy.GetParameterAsText(0)
path_out = arcpy.GetParameterAsText(1)

fname = file_in.rsplit('\\', 1)[1]

###
## Calculate Slope
###

arcpy.env.workspace = path_out
fout = path_out + '\\'+ fname[0:3] + '_slope.tif'
slope = arcpy.Slope_3d(file_in, fout, "DEGREE")
arcpy.env.overwriteOutput = True
files = arcpy.ListFiles('*.tif')

###
## Tiling 2 by 2
###

for i in range(1,4):
    for fi in files:
            fout = str(fi)[:-4] + '_'
            tiles_path = "%s\\%s\\" %(path_in, i)
            if not arcpy.Exists(tiles_path):
                arcpy.CreateFolder_management(path_in, "%s" % i)
            arcpy.SplitRaster_management(fi, tiles_path, fout, "NUMBER_OF_TILES", "TIFF", "BILINEAR", "2 2", "#", "4", "PIXELS", "#", "#")
    arcpy.env.workspace = tiles_path
    files = arcpy.ListFiles('*.tif')

But the following error occured.

  Traceback (most recent call last): File "C:\arcgisserver\directories\arcgissystem\arcgisinput\GeoTools\Slope3.GPServer\extracted\v101\my_toolboxes\Slope_Steepness_v01.py", line 41, in arcpy.CreateFolder_management(path_out, "tiles") File "c:\program files\arcgis\server\arcpy\arcpy\management.py", line 16613, in CreateFolder raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Output Folder: Dataset C:\test\1_proc_data\tiles does not
  exist or is not supported Failed to execute (CreateFolder)

So I gave up on the folder structure and rewrote the script to store everyting into one folder, I already created. Again this worked fine when I was using the Toolbox in ArcMap. But after publishing the Service as WPS, the following error occured:

The second:

import os, sys
import arcpy

####
## Environment Settings ArcPy
####

file_in = arcpy.GetParameterAsText(0)
path_out = arcpy.GetParameterAsText(1)

fname = file_in.rsplit('\\', 1)[1]

###
## Calculate Slope
###

arcpy.env.workspace = path_out
fout = path_out + '\\'+ fname[0:3] + '_slope.tif'
slope = arcpy.Slope_3d(file_in, fout, "DEGREE")
arcpy.env.overwriteOutput = True
files = arcpy.ListFiles('*.tif')

###
## Tiling 2 by 2
###

fout = fname[:-4] + '_'

for i in range(1,4):
    for fi in files:
            fout = str(fi)[:-4] + '_'
            arcpy.SplitRaster_management(fi, path_out, fout, "NUMBER_OF_TILES", "TIFF", "BILINEAR", "2 2", "#", "4", "PIXELS", "#", "#")
    arcpy.env.workspace = path_out
    files = arcpy.ListFiles('*.tif')

  Traceback (most recent call last): File "C:\arcgisserver\directories\arcgissystem\arcgisinput\GeoTools\Slope3.GPServer\extracted\v101\my_toolboxes\Slope_Steepness_v01.py", line 45, in slope = arcpy.Slope_3d(file_in, fout, "DEGREE") File "c:\program files\arcgis\server\arcpy\arcpy\ddd.py", line 4117, in Slope raise e ExecuteError: ERROR 000875: Output raster:
  C:\Users\stankowski\ownCloud\Project_v01\1_proc_data\tiles\srt_slope.tif
  workspace is an invalid output workspace.

I checked the spelling, the use of characters and the length of the output file. I can not make sense of this. The Analysis Tool doesn't find any bugs before publishing it to the server. The only reason I can think of is that there are issues accessing local folders via WPS. Are there any ways to resolve this? I enabled upload to the Server in the Service Definiton File.


texts: Drag and Drop with ArcGIS Engine MapControl within WindowsFormsHost in WPF

texts:
I have an application that requires dragging and dropping onto ArcGIS Engine's MapControl in a Windowsformhost. I have set the Allowdrop property on both the host and the control itself and it still does not allow me to drop on the control.

Is there a way to enable dropping on the MapControl?


texts: Sketch event in Openlayers 3 drawing interaction

texts:
For drawing Polygon in OpenLayers 2 we could use as follow:

function drawEnd(){
    ...
}

function putPoint(){
    ...
}

function sketchModified(){
    ...
}

function cancel(){
    ...
}

var callBackParams = {
    "done": drawEnd,
    "point": putPoint,
    "modify": sketchModified,
    "cancel": cancel
}

var drawCrtl = new OpenLayers.Control.DrawFeature(layer, OpenLayers.Handler.Polygon, {callbacks: callBackParams});
map.addControls([drawCrtl]);
drawCrtl.activate();

When we put point to drawing the polygon, the system call putPoint function. When mouse move over the map, it call sketchModified function. I want to use these two Event in OpenLayers 3. I define my layer as follow:

var source = new ol.source.Vector({wrapX: false});

var vector = new ol.layer.Vector({
  source: source
});
draw = new ol.interaction.Draw({
    source: source,
    type: 'Polygon'
});

drawStart = function(){
    console.log("Salam Bar Mahdi");
}
drawEnd = function(){
    console.log("Ya Mahdi");
}

draw.on('drawstart', drawStart);
draw.on('drawend', drawEnd);
var map = new ol.Map({
    interactions: ol.interaction.defaults().extend([select, modify, draw]),
    target: 'map',
    layers: [
      new ol.layer.Tile({
        title: "Global Imagery",
        source: new ol.source.TileWMS({
          url: 'http://localhost:8080/geoserver/world/wms?service=WMS',
          params: {LAYERS: 'world:worldRaster', VERSION: '1.1.1'}
        })
      }),
      vector
    ],
    view: new ol.View({
      projection: 'EPSG:4326',
      center: [0, 0],
      zoom: 4
    })
});

But in Openlayers 3 I saw only drawEnd and DrawStart Event.

How Can I do?


texts: Changing values in column based on another column in arcpy

texts:
I created new field in a table named 'Cimp'. Now I want to put there values that depends on values from column "COS_Cod" for example, if COS_Cod starts with 1.1 then Cimp is equal to 1. COS_Cod is a text, Cimp must be a number, because I will use it later for calculations.

I write code that base on code that I found in "Refining an Automated Model for Basic Landform Classification A Python and ArcGIS Approach" by Josh Moss. I doesn't give me errors, but when I check table my Cimp is unchanged.

This is my code:

#adding the value
tableFile = "C:/Users/User/Desktop/project/COS_CAOP/results/ClippedCOS.dbf"
arcpy.AddField_management(tableFile, "Cimp", "FLOAT","","","","","NULLABLE")
print ("TABLE DONE")

#Let's try this method!

COS_Cod = "COS_Cod"
Cimp = "Cimp"

with arcpy.da.UpdateCursor(pClippedCOS, (COS_Cod, Cimp)) as updateRows:
    for lForms in updateRows:
        cosCodes = lForms[0]
        print cosCodes

        print "printed cosCodes"

        def codeCimpAssign(COS_Cod):
            if cosCodes == "1.1*":
                numCimp = 1.0
            elif cosCodes == "1.2*":
                numCimp = 1.0
            elif cosCodes == "1.3*":
                numCimp = 0.8
            elif cosCodes == "1.4*":
                numCimp = 0.4
            elif cosCodes == "2.1*" and cosCodes == "2.2*":
                numCimp = 0.35
            elif cosCodes == "2.4*":
                numCimp = 0.35
            elif cosCodes == "3.1*":
                numCimp = 0.0
            elif cosCodes == "3.2*":
                numCimp = 0.25
            elif cosCodes == "3.3*":
                numCimp = 0.8
            else:
                numCimp = 0.0
                print "something"
            return numCimp

# define "text" as relationship for updating rows
        numCimpValue = codeCimpAssign(cosCodes)
        print numCimpValue
        print "printed Cimp values"

#update rows

        lForms[1] = numCimpValue
        updateRows.updateRow(lForms)
        print "Added CImp values"
        arcpy.AddMessage("Added CImp values")


texts:
I have 4 vertices and I need to get the bounding box something like the one in the image below.
A, B, C and D are the vertices.

The red area is the one which I need to eliminate while I only need to retain the green area.

Maybe, I would need ST_Intersects, but I'm not sure how do I go about using this for the mentioned use.

NOTE I will need the points A, B, C and D as part of the green area, I do not want to eliminate them, I only need to get rid of the center part which is in red.


texts: How to create a color ramp in ArcGIS 10 with only 10 colors

texts: I am trying to use the company I work for's color scheme in the maps that I am making.  I have all the RGB values for the colors, but there are only 10 different shades of each color, while the style manager's ramps need 12 shades.  Is there any way to change the number of colors in each color ramp?


texts:
On my way to find a solution to This Question, I´m trying to merge some automatic generated script with This example.

My thought is that after the last step, I could insert the script that does nothing with the file but only rename the output temporary file, because it´s not a native source.

Running everything in the toolbox.

How can I set the temporary output as an input to the example script?

This is the part of the automatic generated script:

        OUTPUT='OUTPUT'

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('maininput', 'Main Input', types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT,self.tr('My output layer')))

        # Buffer
        alg_params = {
            'DISSOLVE': False,
            'DISTANCE': 10,
            'END_CAP_STYLE': 0,
            'INPUT': parameters['maininput'],
            'JOIN_STYLE': 0,
            'MITER_LIMIT': 2,
            'SEGMENTS': 5,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Buffer'] = processing.run('native:buffer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

and after that, I want to insert the other script, having "outputs['Buffer']['OUTPUT']" instead of "INPUT"

        source = self.parameterAsSource(parameters, outputs['Buffer']['OUTPUT'], context)

        (sink, dest_id) = self.parameterAsSink(
            parameters, "OUTPUT", context,
        source.fields(), source.wkbType(), source.sourceCrs())

        total = 100.0 / source.featureCount() if source.featureCount() else 0
        features = source.getFeatures()
        for current, feature in enumerate(features):
            if feedback.isCanceled():
                break
            out_feature = QgsFeature(feature)
            sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
            feedback.setProgress(int(current * total))

        return {"OUTPUT": dest_id}

I got an error

  AttributeError: 'NoneType' object has no attribute 'fields'


texts:
I'll prefix this with stating that my abilities in GIS are lackluster and I'm still learning, it is also my first time doing any 3D work in QGIS.

I have a KMZ file that contains Z Values (Created in Autocad) which I would like to use in a 3D GIS file.

Will this require use of the "Mesh Calculator"?


texts: I can't find TauDEM in the providers under processing

texts: I am trying to add TauDEM, after installing it as instructed I can't find it Under the providers in QGIS 3.4 although I can under earlier versions (also can't install them for other reasons)


texts: Generating perfect tile set from Mapbox for given zoom level?

texts:
I'm using Mapbox API to generate bitmap tiles.
First of all I have generated whole world tile with zoom 1 (seems it doesn't matter 0 or 1) and resolution 1280x1280. This tile perfectly positioned from -180 to 180 lon.

Next I would want to generate 4 tiles for next zoom 2. If I correct first tile would be at -90.0,90.0 position (around North America).

But I noticed that size of resulted tile is not as I wanted. It is bigger.

So how to found appropriate zoom level or something else to get tiles for my grid like this, and then for factor 4, 8 etc.


texts: Open popup dialog window on button click in Python AddIn?

texts:
I'm trying to figure out a way to display an image in a dialog box that will be used for a help menu for my toolbar addins.

The best scenario i've found is using the pythonaddins module and the messagebox popup, but can't figure out how to any formatting, except pure text string.

Is there another way to display information that i can call via a button click.


texts:
I have a CARTO map displayed via Leaflet (specifically via Leaflet within the CARTO.js library) using carto.createlayer().

I simply want to get the_geom values from each marker as they're clicked on. Layer interactivity is set to the_geom but values are being returned as null. I can get the code below to return the values of other columns in my dataset (such as cartodb_id, site_id, etc.), but I don't seem to be able to access the_geom values using data.the_geom.

cartodb.createLayer(map, {
  user_name: 'codemacabre',
  type: 'cartodb',
  sublayers: [{
    sql: 'SELECT * FROM hydrous_sites_geo_map',
    cartocss: '#hydrous_sites_geo_map {marker-fill: #F0F0F0;}',
    interactivity: 'cartodb_id, the_geom'
  }]
})
.addTo(map)
.on('done', function(layer) {
  layer.setInteraction(true);
  layer.on('featureClick', function(e, latlng, pos, data) {
    console.log(data.cartodb_id);
    console.log(data.the_geom); // returns null
  });
});

See the screenshots for the console output (top) and actual dataset values (bottom):

I've read that the_geom can become corrupted when exporting / re-importing (such as in these issues: 1 & 2), but I'm unsure if that's the case here or if the_geom protected somehow or (most likely) if I'm missing something obvious.

Is there a knack to returning the_geom values?


texts: Communicating with Java Class from OpenLayers 3 UI?

texts:
I have an OpenLayers 3 map and I am using Geoserver to host that map. So far, I have mostly frontend related tasks for my mapping application. Now, I need to include zonal statistics feature in my application. I was looking online and found that GeoServer/GeoTools RasterZonalStatistics.Java file is what I need to compute my zonal statistics.

However, I don't know where to write my zonal statistics class.

How can I communicate with this Java class from OpenLayers UI?

For example, should I write my zonal statistics class in a separate eclipse project and host that Java file into GeoServer?

I tried to look online and came across GeoServer's OWS services. However, I am not sure if that's what I need to set up a correct backend class for my OpenLayers app.

How can I communicate with this Java class from my OpenLayers UI to calculate zonal statistics for my Tiff and shp file?


texts: An error occurred when using the terrain function in the production of Shiny dashboard

texts:
I am new to the 'shiny' package. I just try to produce an aspect map from the existing DEM data using the Shiny dashboard for training purposes. However, I faced the below warning in R and error in Shiny dashboard.
Warning: Error in h: error in evaluating the argument 'x' in selecting a method for function 'plot': unable to find an inherited method for function ‘terrain’ for signature ‘"NULL"’
Error:error in evaluating the argument 'x' in selecting a method for function 'plot': unable to find an inherited method for function ‘terrain’ for signature ‘"character"’
I realized situations where similar errors were encountered but I didn't find the exact reason and didn't solve this issue.
Similar error - I: What does this mean: unable to find an inherited method for function ‘A’ for signature ‘“B”’
Similar error - II: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘shapefile’ for signature ‘"NULL"’
Here is my simple code:
library(shiny)
library(raster)
library(rgdal)

ui=shinyUI(fluidPage(pageWithSidebar(
  headerPanel("Header1"),
  sidebarPanel(
    fileInput('DEM', 'Choose DEM', multiple=FALSE, accept='asc'),

  ),
  mainPanel(
    plotOutput("aspect")
  )
)))

server = shinyServer(function(input,output){

  inFile <- reactive({
    raster::terrain(input$DEM$datapath, opt = 'aspect')
  })

  output$aspect<-renderPlot(
    {
      plot(inFile())
    })
})

shinyApp(ui = ui, server = server)


texts:
The problem consists in modelling propagation of a signal (eg. light or sound etc) through a series of obstacles, as on the figure below. The signal cannot pass through bottom surface (terrain), but it can traverse obstacles. I want to count the number of obstacles traversed.

Terrain and obstacles are in 2D numpy arrays (x,y,z). This is what I do:

output = numpy.zeros(terrain.shape)

obstacles = terrain + obstacle_heights

for i in xrange (obstacles.shape[0]):
    for j in xrange (obstacles.shape[1]):

        mask = obstacles[i,j] > terrain[i,j:]
        output[i,j:][mask] +=1

The result would be something like [0, 0, 0, 1, 1, 1, 2, 3, 4, 4, 4 ...] per row.

This method works fine (provided valleys on the terrain are filled by using numpy.maximum.accumulate). Now, would it be possible to speed up the thing by using a vectorised solution?


texts: Add new attribute in shapefile using Java: Geotools

texts:
I want to edit my shapefile by adding a new column "ID" and assign a value to it. I am trying to achieve this using Java Geotools. Here is the snippet code:

 File inFile = new File("/abc".shp");

        System.out.println("********* GEOTOOLS *************");

        System.out.println("********* READ *************");
        Map<String, Object> map = new HashMap<>();
        map.put("url", inFile.toURI().toURL());

        DataStore dataStore = DataStoreFinder.getDataStore(map);
        String typeName = dataStore.getTypeNames()[0];

        FeatureSource<SimpleFeatureType, SimpleFeature> source =
            dataStore.getFeatureSource(typeName);
        Filter filter = Filter.INCLUDE; // ECQL.toFilter("BBOX(THE_GEOM, 10,20,30,40)")

    FeatureCollection<SimpleFeatureType, SimpleFeature> collection = source.getFeatures(filter);
    try (FeatureIterator<SimpleFeature> features = collection.features()) {
        while (features.hasNext()) {
            SimpleFeature feature = features.next();

            // Adding new column attribute ID

            feature.setAttribute("id", "ADSAS2132312");

            System.out.print(feature.getID());
            System.out.print(": ");
            System.out.println(feature.getDefaultGeometryProperty().getValue());
        }
    }

// Write the features to the shapefile

System.out.println("********* WRITE *************");

Transaction transaction = new DefaultTransaction("create");

   // String typeName = newDataStore.getTypeNames()[0];
    SimpleFeatureSource featureSource = newDataStore.getFeatureSource(typeName);

    if (featureSource instanceof SimpleFeatureStore) {
        SimpleFeatureStore featureStore = (SimpleFeatureStore) featureSource;

        featureStore.setTransaction(transaction);
        try {
            featureStore.addFeatures(collection);
            transaction.commit();

        } catch (Exception problem) {
            problem.printStackTrace();
            transaction.rollback();

        } finally {
            transaction.close();
        }
        System.exit(0); // success!
    } else {
        System.out.println(typeName + " does not support read/write access");
        System.exit(1);
    }
}

I am not able to write the ID column into shapefile. The reading works fine, but not the write part. I don't know much of Java, so its hard for me to figure out. I followed this reference:
http://docs.geotools.org/latest/userguide/tutorial/feature/csv2shp.html


texts: Add and remove tools from ArcMap toolbar via C# Button and ArcObjects?

texts:
I am trying to create a custom button in arcmap that when pressed will edit a toolbars content. Doing things like adding or removing tools from a toolbar. An example would be to create a button that when pressed removes or toggles all the text/font editing tools from the draw toolbar.

I've been Looking at with IToolbarControl2
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#/IToolbarControl2_Interface/0016000003m7000000/
It seems to give access to various member functions that let me edit what I want. unfortunately I'm unsure of how to reference the specific toolbar that I want to edit.

I made some edits and currently the code runs and then hangs when i attempt to define axToolbarControl1.AddToolbarDef via:

axToolbarControl1.AddToolbarDef(uID, -1, false, 0, esriCommandStyles.esriCommandStyleIconOnly);

The current draft of the code is the following:

///Code to add or remove tools from a toolbar using System;
using System;
using System.Collections.Generic;
using System.Text;
using System.Windows.Forms;
using System.IO;
using ESRI.ArcGIS.Framework;
using ESRI.ArcGIS.ArcMap;
using ESRI.ArcGIS.ArcMapUI;
using ESRI.ArcGIS.Controls;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.SystemUI;
using ESRI.ArcGIS;

namespace RemoveAllTools2
{
public class RemoveAllTools2 : ESRI.ArcGIS.Desktop.AddIns.Button
{
    IApplication m_application;        //ArcMap application
    IMxDocument m_mxDocument;          //ArcMap document
    public ESRI.ArcGIS.Controls.ToolbarControl axToolbarControl1;
    public UID uID = new UIDClass();

    public RemoveAllTools2()
    {
        MessageBox.Show("setup start");
        m_application = ArcMap.Application;
        MessageBox.Show("m_application set");
        UID uID = new UIDClass();
        uID.Value = "esriArcMapUI.DrawToolBar";
        MessageBox.Show("UID set");

        //code hangs here
        axToolbarControl1.AddToolbarDef(uID, -1, false, 0, esriCommandStyles.esriCommandStyleIconOnly);
        //code hangs here

        MessageBox.Show("axtoolbar set");

    }

    protected override void OnClick()
    {
        MessageBox.Show("clicked");
        IToolbarControl2 ttoolbarControl = axToolbarControl1.Object as IToolbarControl2;
        MessageBox.Show("toolbarControl set");
        uID.Value = "esriArcMapUI.FontNameControl";

        int toolFindIndex = ttoolbarControl.Find(uID);
        if (toolFindIndex!=-1)
        {
             ttoolbarControl.Remove(toolFindIndex);
             MessageBox.Show("Tool remove run finished ");
        }
        else
        {
            MessageBox.Show("Tool not found");
        }
        MessageBox.Show("IF run finished ");

    }
    protected override void OnUpdate()
    {

    }
}

}


texts: Identify Health Centre Location using machine learning ( Optimisation )

texts:
I am trying to identify locations of where new health centres could be built in the country using Machine Learning (optimisation) and R programming language. I currently have the latitude and longitude of current health centres. I was looking at this tutorial:

https://www.dataplusscience.com/GeoOptTabInstructions.html

and I'm not sure if its adequate for what Im trying to do. I was also looking at the Bayesian algorithm as well.

Any thoughts on which way to proceed?