texts: How to use gdal_grid with ungridded .asc DTM files

texts:
There was a question about processing ungridded .asc DTM files: Problem processing Slovenian LiDAR DTM with GDAL (and QGIS): "too many stepY values" error
The answer helped a lot already, but I don't understand how you calculate -txe and -txy for the gdal_grid command?
When I look into the example file "GK1_444_106.asc", the first point (lower left) is:
444000.00;106000.15;799.16
and the last point (upper right) is:
444999.96;106999.15;712.21
I would assume the georeferenced extents of the gdal_grid command would be:
-txe <xmin> <xmax> => -txe 444000.00 444999.96

-tye <ymin> <ymax> => -tye 106000.15 106999.15

But in the answer it is:
-txe 443999.5 445000.5 -tye 106999.65 105999.65

Why was there 0.5 added (to max) and subtracted (from min) to the values (and why were they rounded)?
Why are the -tye values in a reversed order (max first, then min)?


texts: Passing attributes from one shapefile to another that does not share line or segment?

texts:
I am struggling to convert attributes from one shapefile to another. The issue is that the geometry is not the same.

How is it possible to pass attribute values from one layer to another considering they don't match ? See below an example.


texts: Stretching a map to the full width of the browser after complete zoom out and without repetition of grids?

texts:
I want to extend the first block from left (with red inside borders) to stretch to the full width of the website. The example code is given below.

html, body {
  height: 100%;
  background: #E6E6E6;
  margin: 0;
  font: 10px sans-serif;
}
svg {
  display: block;
}
circle {
  stroke: black;
  fill: brown;
  fill-opacity: .5;
}
#map {
  width: 1370px;
  height: 600px;
}
    </style>


texts:
I try to add a WMS layer to QGIS Desktop from my plugin (C++)

QString uri = "crs=EPSG:28992&layers=natura2000:natura2000&styles=&format=image/png&url=http://geodata.nationaalgeoregister.nl/natura2000/ows?
QString myProviderName = "wms";
QgsRasterLayer* lyr = mQGisIface->addRasterLayer(uri ,"test", myProviderName);

But in QGIS I get the message, the extent of the layer could not be calculated because of unknown dataprovider.

What am I doing wrong?

In QGIS Desktop manually adding this WMS layer is possible without problems.
I also tested several other WMS with same result.


texts: What is the difference between QueryTask and featureLayer.queryFeatures

texts:
For fetching data from a feature layer, I've been using QueryTask service,

    var deferred = new Deferred();
    var queryTask = new QueryTask(featureLayerurl)
    var query = new Query();
     query.geometry = geometry;
    query.where = "document='100'";
    query.returnGeometry = true;
    query.outFields = ["*"];

    return deferred = queryTask.execute(query).addErrback(function (response) {
        console.error(response);
        return null;
    }).addCallback(function (response) {
        return response;
    }).then(function (response) {
        return response;
    });

But recently, I've found that there is dedicated method in FeatureLayer class called queryFeatures for for fetching data.

var query = new esri.tasks.Query();
query.objectIds = [features[0].attributes.OBJECTID];
query.outFields = [ "*" ];
// Query for the features with the given object ID
featureLayer.queryFeatures(query, function(featureSet) {
});

So ,

Is there any difference between when using QueryTask over
queryFeatures?
Any  kind a of performance difference ?
Do queryFeautres is using QueryTask internally ?


texts:
I would like to open a Plugin with the help of another Plugin.
Is it possible, if I want to switch the workflow of different plugins one after the other, that when I click the submit button of one plugin I close this plugin and another one opens automatically?


texts:
I am already using proj4j in my project to convert between different coordinate reference systems. (Internally I use WGS84 but some sources supply data in a different CRS.)

Now I need to calculate the distance between two points in meters and am wondering if (and how) that is possible with PROJ4J. (Just in case: my use case is limited to a maximum distance of around 100 km.)

I have come across org.osgeo.proj4j.util.ProjectionMath#greatCircleDistance(double, double, double, double), but that class is not documented at all. What CRS and unit does that method expect for the latitude/longitude values? And what is the unit for the result?


texts: Get band number of maximum value per pixel for multiband raster in R

texts:
I have a multiband raster and I want to find the band number of the maximum value for each pixel.

Here is some example code:

library(raster)
set.seed(42)
r1 <- raster(nrow=4, ncol=4)
r1[] <- round(runif(ncell(r1)) * 10)
r2 = raster(nrow=4,ncol=4)
r2[] = round(runif(ncell(r2)) * 10)
values(r1)
values(r2)

rast = brick(r1,r2)

Obviously, I can find the maximum value for each pixel, but I want to know from which band it came. In the above example values(r1) = [1] 9 9 3 8 6 5 7 1 7 7 5 7 9 3 5 9 and values(r2) = [1] 10  1  5  6  9  1 10  9  1  5  4  9  4  8  7  8

max(rast) = [1] 10  9  5  8  9  5 10  9  7  7  5  9  9  8  7  9

So I want a raster that eventually looks like [2 1 1 1(2?) . . . . 1].

How do I do this?


texts: Automatically adapt label text size to area in QGIS

texts:
In Relating text buffer to text size in QGIS? I displayed various polygon sizes : the objective of my research was to try to adjust automatically the buffer size of the text to the font size because the buffer tended to be too large for little texts.

I'm now looking for a way to adapt the font size to the polygon area. This way, i won't have to put it into the attribute table and fix manually the font size i need.


texts: ArcMap python add-in button fails to load and displays Missing on the toolbar

texts:
I have a python add-in button shown as follows. It has to make a http request and use the values from the response. But whenever I use requests package, the add-in result in Missing button.

import arcpy
import pythonaddins
import requests

class ButtonClass1(object):
    """Implementation for ZoomButton_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        url = "https://www.linearbench.com/api/azgiv/extent"
        r = requests.get(url = url, verify = False)
        #data = r.json()
        #xmin = data['xmin']
        #ymin = data['ymin']
        #xmax = data['xmax']
        #ymax = data['ymax']
        mxd = arcpy.mapping.MapDocument('current')
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        #ext = arcpy.Extent(xmin,ymin,xmax,ymax)
        ext = arcpy.Extent(45,2,50,5)
        df.extent = ext

According to Esri, it's commonly caused by syntax error.

  Anytime a tool or button fails to load and displays Missing on the toolbar instead of the expected icon or name, it is commonly due to a Python syntax error. Recheck your Python script for syntax errors or incorrect inputs to Python functions.

But I tested my code in the ArcMap python window. It works fine. What could be the cause of the problem?


texts: GDAL 2.4.2 installation problem in Debian

texts:
I have built GDAL 2.4.2 in Debian package through the following commands

sudo apt-get install build-essential python-all-dev
wget https://download.osgeo.org/gdal/2.4.2/gdal-2.4.2.tar.gz
tar xvfz gdal-2.4.2.tar.gz
cd gdal-2.4.2
./configure --with-python
make clean
make
sudo make install
sudo ldconfig

No error shows during building. Now ogr2ogr command is not working and returning -bash: /usr/bin/ogr2ogr: No such file or directory. Whereas it is showing that gdal is installed if I asking for the version gdalinfo --version. It's returning GDAL 2.4.2, released 2019/06/28.

How can I get rid of this in my Debian server?


texts: Real-world dimensions of a 3D object in VisualSFM or MeshLab

texts:
I have a good-quality, precise 3D scene reconstructed from photographs/video by VisualSFM. All camera parameters, including 3D position and lens details, are known.

Now, that the hardest work has been accomplished, I'm sure that obtaining real-world dimensions of the model is easy... only how?

Nobody on the Internet besides me appears to want to know – no answers to this elementary problem exist on the web.

If VisualSFM lack necessary function, I can use an external tool. Example: MeshLab.

If no extant software has necessary function, it shouldn't be too hard to write a C programme or an Excel spreadsheet formula to this end... if someone can tell me all needed formulas.


texts: ArcGIS REST API find features which contain given point

texts:
I've developed an application which loads the services from the ArcGIS Server via REST API.

I load a service, the service has multiple layers, each layer has multiple features (objects).

By the URL http://domain/arcgis/rest/services/Test/MapServer/3/query?f=json&where=1%3D1&returnGeometry=false&outFields=* I get all features of a layer.

When I click on the map, I get the Lat/long of the point. I want to load all features which contain this point. How should I generate the URL?


texts: QGIS Labeling OSM items depending on scale

texts:
I'm working with OSM data in QGIS 2.10 (Ubuntu). I'd like to label items depending on scale. What I mean to do is to have

only labels of "primary" roads from a scale of 1:80.000 down to 1:50:000
all the other labels from 1:50.000 down to 1:1.000

The only expression about labeling that sounded helpful was:

CASE WHEN "place" in ( 'city','town' , 'village' , 'suburb' )
THEN 0
ELSE 1
END

But it doesn't exactly do what I'm looking for.

Does anybody know what I should do?


texts: Executing standalone Python script for QGIS processing from within PowerShell script

texts:
Is there a way to execute either a python script or pyhton commands from within a powershell script to do some QGIS geoprocessing on a Windows 10 machine (64 bit)?
I have a powershell script that downloads gazetteer data with cURL and copies the data into a geopackage. I want to enrich the address data with a number of other attributes from other layers. As of now, I use python processing commands for this task that need to be executed from the python console of QGIS. A third script (powershell again) exports the data to csv files, 7z-ips and copies them to the clipboard in preparation for mail delivery.
Example of processing:
processing.run("native:joinattributesbylocation", {'INPUT':'C:/Users/myusername/Documents/Parkrauminventur.gpkg|layername=Hauskoordinaten','JOIN':'//some_server/Projekte/PRM.gdb|layername=Kontrollgebiete','PREDICATE':[0,1,2,3,4,5,6],'JOIN_FIELDS':['Gebiet','Region'],'METHOD':1,'DISCARD_NONMATCHING':False,'PREFIX':'KG_','OUTPUT':'ogr:dbname=\'C:/Users/myusername/Spatial_Joins.gpkg\' table=\"SJ_HK\" (geom)'})

However, I would like to write one script to include all steps that runs without the need to start QGIS. I read many instructions on how to write standalone scripts but none of them seem to be assuming powershell. So I am not sure how to set the required environmental variables in powershell and how to interact with python/QGIS.
Are there directions for this scenario?


texts: Export NDVI values from chart as JSON in Google Earth Engine

texts:
After filter OLI8 for a small area

var collectionFilter = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA')
  .filterBounds(geometry)
  .filterMetadata('CLOUD_COVER','less_than',4)
  .map(app);
print(collectionFilter)

And add NDVI band

var app = function (image)
{
//indice de vegetação
var ndvi = image.normalizedDifference(['B5', 'B4']);
image=image.addBands (ndvi.rename('NDVI'));

  return image
}

I generate a chart using a click function and try to export as JSON format the values.

Map.onClick(function(lonlat){
  var pixel = ee.Geometry.Point([lonlat.lon, lonlat.lat])
  var collection = collectionFilter
    .select(
      ["NDVI"]
    )
    var options = {
        title : 'NDVI CHART VALUES'
      };

    var graph = (
    ui.Chart.image.series(
      collection,
      pixel,
      ee.Reducer.mean(),
      30
    )
    .setOptions(options)
  );
      print(graph

  )

  Export.table.toDrive({
    collection: collection,
    description:  'exporting in json format each click',
    fileFormat: 'JSON'
  });

})

But doesnt work. How can I do this? Here is the link code


texts: How to route between lat/lon pairs in pgRouting?

texts:
I have set up Postgresql, PostGIS and pgrouting on Ubuntu.
I have imported an osm for wales using the following command:

./osm2pgrouting -file wales.osm \
                            -conf mapconfig.xml \
                            -dbname routing \
                            -user postgres \
                            -clean

All I now want to do is execute something to give me the fastest route between 2 Latitude/Longitude points.

Start: 51.499954,-3.571651
End: 53.386809,-4.445

The route has to be in the format of a series of latitude/longitude positions that make up the route.

Is this even possible? If not, what would it take to make it possible.

If it is possible, what sql do I need to execute?

I have been looking everywhere for an answer to no avail.

Further to Steve Horn's recommendation. I have now managed to get 5000 points coming back using this sql:

SELECT gid, AsText(the_geom) AS the_geom FROM dijkstra_sp('ways', 15522, 26857);

However, the points all seem to be ignoring the road type.
(When I run it using Google Directions, the route mainly follows the main roads).
What have I done wrong?


texts: Why are size-scaled symbols missing from map layout legend?

texts:
I have a designed a map using a scaled simple marker symbol. It is a rule-based design and symbol scaling is done by the siye assistant based on a field value for any of the rules in the same way. Everything works fine within the the map - legend / table of contents symbols are scaled according to map scale, within the print composer / item properties-window symbol scaling works fine as well.

However, the legend of my map-layout doesn't show any symbols.

Strange: changing scales in the layout (via  map/item properties/scale) changes line spacing of my map-legend - indicating that there is information about how to scale symbols somewhere in the background.

... but the symbols remain invisible in the map-legend. Does anyone have any ideas?

QGIS 3.2 Mac and Windows.

@tallistroan
seems that there is a serious problem in version 3.2 ...
following your workaround description resulted in correct scaled symbols on the map and ugly results in both the Advanced > Data-defined Size Legend tool and the composer legend. obviously the sizes do not correspond to the map scale and to symbol sizes within the map.

qgis application window / advanced data-defined size tool:composer output:

composer window - ledend item properties and map-output:


texts:
I have a XML file that would like to use as a base layer in ArcGIS 10. What are the steps that I need to follow in order to open that file in ArcGIS 10?

The XML file does not have coordinates. When it's open it shows the following:


texts:
I hope somebody has already figured this one out. I just installed Geoserver 2.9 on a vanilla Ubuntu 16.04 distro. The Geoserver 2.8 method of enabling CORS with the shanbe.hezoun class does no longer work with Jetty 9.2.13.

There are mentions that CORS support is already packaged with Jetty 9.2.13 in the jetty-servlets.jar.

The Jetty lib which is compiled with Geoserver contains a jetty-servlet-9.2.13.v20150730.jar in geoserver/lib but not jetty-servlets.9.2.13.v20150730.jar. Are these supposed to be the same jar with a different name?

It should be possible to enable CORS either in geoserver/etc/webdefault.xml or in geoserver/webapps/geoserver/WEB-INF/web.xml.

My understanding is that the webdefault.xml is applied first and the web.xml thereafter.

I have tried following filter in both xml. I haven't got as far as adding a filter mapping. Adding the filter alone will cause the Geoserver/Jetty service to not start proper.

<filter>
    <filter-name>cross-origin</filter-name>
    <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
</filter>


texts: Adding field values from feature class into empty field of another feature class using append

texts:
I am attempting to add a feature class ("dis.shp") into an empty feature class("fc.shp") in a geodatabase using Python script.  I have had success with the code below.  However, I also would like to add values from dis.shp's field "GRIDCODE" into fc.shp's field "Classes".  I have been playing around with append but am not having any luck.  Any suggestions?

arcpy.Append_management(dis, fc , schemaType, fieldMappings, subtype)


texts: Calculate lengths by mode in network analyst

texts: I have a multi-modal network (Rail, Roads and River/Ocean) for network analyst.  I am calculating travel costs along the multi-modal network.  I get a total length field, but is there anyway to get a total length by each mode for each route?  I have thousands or routes so don't really want to try selecting them all out.  Any ideas?


texts: How to use ModelBuilder to add many images to mosaic dataset?

texts:
Here is my sample of a model builder. What I would like to create a mosaic dataset then add lots lots lots lots of imagery from different folders to the dataset rather than I have to do manually.

What I would like to run tool to add the imagery to the dataset and when it is done then change to another folder and keep adding more imagery without having to stop the process builder. Is there a tool to do that or should I need to create a variable ?

There are three folders that hold imagery and there are three sub-folders as well.

You can see why I am trying to build three different names: Columbus, El Paso, and Playas Valley because it is where they hold the imagery.


texts: How add a text field and then reproject all shapefiles which are in sub folders

texts:
I have a folder and inside of that folder, I have about 500 sub folders which contain on average about 5 shapefiles each.

What I want to do is:

to check if there is a field, lets say field_1, and if it is not in  the attribute table, to add it (TEXT length = 6) and then
to reproject all the shapefiles preferably with ogr2org to EPSG:3035

I find this question 1 but my data are in multiple directories and this question 2 which is for converting all shapefiles in a directory.

the above questions give me food for thought of how can someone do something that I am looking for.


texts: Drawing geometrical shapes with dimensions using ArcGIS API for JavaScript?

texts:
I have to draw the geometrical shapes with their dimensions and show it with that shape.

As in below images,

I want to show the length and width of the rectangle, the radius of the circle and the area of the shape at centre.
If the user selects and edits the geometrical shapes then the dimensions should be updated accordingly.

I have read the below examples, but I couldn't able to find the proper one.

Add graphics to a map | ArcGIS API for JavaScript 3.25
Drawing tools | ArcGIS API for JavaScript 3.25
Measurement Widget | ArcGIS API for JavaScript 3.25
Reshape polygons | ArcGIS API for JavaScript 3.25


texts:
I am trying out Mapbox-gl-js as an alternative to our current openLayers 3 solution. This is primarily due to a desire to shift away from raster backgrounds and to utilize the webGL and vector tile support in Mapbox-gl-js. One issue i've encountered and am not exactly sure how to tackle is how to maintain and accurate and scaling vector curve in Mapbox-gl-js (perhaps using turf.js?). I have not been able to find an example else where so hopefully someone out there brighter than me can help. In the example below you will see the difference between Openlayers3 and Mapbox-gl-js:

OpenLayers 3:

Mapbox-gl-js:

Currently the data is brought into the map via an external geoJSON file and just added to the map with a simple call like so:

     //Define source layer
     map.addSource('lineSource', {
        type: "geojson",
        data: lineGeoJSON //geojson path
     });

     //Define Mapbox layer for line
     map.addLayer({
        "id": "lineGeoJSON",
        "type": "line",
        "source": "lineSource"
     });

FIX: Not 100% sure of how I fixed this. I switched from an ajax to an xml call to grab and parse my geojson data. I also began converting my web mercator meter coordinates to lat/long on the fly in the geojson parser vs. in QGIS. This seems like a rare and odd bug, but hopefully this info may help someone down the line.


texts: How to sum two fields in Python Advanced Field Calculator?

texts:
I need to perform a quite simple task (or at least it was simple on ArcGIS):
I have two numeric fields called c0g1 and c0g2, and want to add them into a new field named TOTAL.
Must use the Advanced calculator because I need to batch the process. I simply do not know what to type in the Formula field (value =). Tried many things but nothing works...


texts: Combining and extracting attributes from 2 columns to populate a third (Concatenate)

texts: I am working with a FEMA flood zone shapefile which contains a FIELD_ZONE field and a SUBTY field.  I created a third field called ZONE_TYPE that needs to be populated by combining or extracting elements of the other two. So if FLD_ZONE = A, Ae, D, V, etc, then ZONE_TYPE needs to be 1 PCT;  if the FIELD_ZONE = X and SUBTY = Annual .2, then ZONE_TYPE needs to be .2 PCT.  I'm guessing I use the field calculator, but I am not a code person (it's like a foreign language to me), and have no idea how to build the block.


texts: How can I load strictly tiles within the current viewport in Leaflet.js?

texts: I have tiles that are not available for the whole map.  How can I stop Leaflet.js/Mapbox.js from making requests to these tiles?  There are dozens of HTTP 404's that I'd really like to get rid of.


texts:
I am trying to create a customized tool in ArcGIS using Python Toolbox. Does anyone know how I can change the appearance (for example, size of a textbox or position of a parameter) of a parameter on the tool dialog? I did some searching online and could only found this ESRI documentation: http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/documenting_tools_and_toolboxes/changing_tool_dialog_appearance.htm, which does not provide much information.
In the built-in ArcGIS Spatial Analyst Tools->Interpolation->Kriging tool, the dialog looks like this:

How can I make my customized tool to have something like the red part in the Kriging tool dialog?
I am using ArcGIS Desktop 10.7.


texts: Style edges of layer shapefiles to have thicker borders

texts:
I have two layers I'm working with here, both of North Carolina. One is "townships" the other is "counties". I want the counties shapes to have thick edges with no fill-in color.

I couldn't find anything on this after a google search so decided to ask.


texts: Linking Coded Domain to Featureclass Field in Enterprise Architect?

texts:
Im using Sparx Enterprise Architect (Ver. 13) to create an ArcGIS Geodatabase.

I have defined a Coded Value Domain but how do I link this Domain to a field within my Featureclass?

I watched the intro video, searched the docs.


texts: Combining two shapefiles into third shapefile in ArcGIS

texts: Suppose I have two shapefiles (A and B). Each contain 200 polygons. There are 150 polygons that are common i.e same polygons in both shapefiles but 50 polygons in each shapefile are unique. I want combine both these shapefiles into third shapefile (C) containing 250 polygons. 150 (common) + 50 (shapefile A) + 50 (shapefile B). I have tried the merge tool but the output shapefile has less number of polygons.


texts: CSWcataloguesource.js plugin "add to map" button not working!

texts:
I got stuck in using CSWcataloguesource plugin in my OpenGeoSuite custom application.
It shows all the layers available in geonetwork but when click the "add to map" button it shows the following error.

TypeError: c.getLayer is not a function

and m not able to add it as a overlay.
Does anyone have idea what this error?and does anyone have used this CSWcataloguesource plugin?

Thanks in advance.


texts: Transposing data points to line using ArcGIS Desktop

texts:
See sketch below.

The main objective, using ArcGIS 10 or 10.1, is to get the P points to the perpendicular P line and, either the closest B point to the endpoint of the P line where it interesects the B lines, or an averaged value, whichever is easier.  My process for doing this is:

"Step 1" is the data as I initially have it.  The parallel B lines have data points and the perpendicular P line has data points scattered around it.

To get to Step 2, I am creating perpendicular lines along the P line then joining the P data point attributes to the perpendicular lines.  Then, I am creating points at the intersection of the P line and the perpendicular lines.  I am then joining the perpendicular line attributes to the intersection points, resulting in the "Step 3" figure.

I am then creating points like those shown in Step 4 and doing a spatial join that assigns the closest value from the points in Step 3 to the F data points in step 4.

Is there a more efficient way of accomplishing this process?

I have about 400 P lines and anywhere between 200-1000 P points at each P line.  I have tried using the Geospatial Modeling Environment, previously Hawth's Tools, with the snap points to line tool but am unable to get it to work "Command text could not be interpreted error).  I think this is because I have P lines that are close together (the sketch below is a simplified idea of what I am doing) and when the tool is trying to snap points to line, the point lies within the snapping tolerance distance of 2 P lines, causing it to crash.


texts: SQL query works in PostgreSQL but not in GeoServer SQL View

texts:
I can't find a solution for this problem after a couple of hours of research. The query is:

SELECT searchaddress('some string');
FETCH ALL IN mycursor;

The function searchaddress returns a REFCURSOR, named mycursor. In PostgreSQL using FETCH I get results. Using the exact same query (no parameters) in a Geoserver SQL View I get an error:

ERROR: syntax error at or near ";" Position: 73

If you need more information I will provide.

Edit: One more piece of information. The query used to be:

SELECT searchaddress('some string');
FETCH ALL IN "<unnamed cursor 1>";

Under this format it was anyone's guess what was the number of "unnamed cursor i". But the SQL View did return a response and not an error when I managed to guess the number of the cursor. Because I don't know how to guess the number I had the function searchaddress changed so it would return a cursor named "mycursor". At that point the querry stoped working altogether under Geosever.


texts: Define Shapefile Projection with Shapely/Fiona

texts:
I am writing a fairly straight forward script of converting GeoJSON data into a shapefile using Shapely and Fiona. My script works but I cannot figure out how to define the CS in the script.

def GeoJSONToFC(shpname,USGSurl):
    import urllib2
    import json
    from shapely.geometry import Point, mapping
    from fiona import collection
    from fiona.crs import from_epsg
    crs = from_epsg(4326)
    schema = {'geometry': 'Point', 'properties': { 'Place': 'str', 'Magnitude': 'str' }}
    with collection(shpname, "w", crs, "ESRI Shapefile", schema) as output:
        url = USGSurl
        weburl = urllib2.urlopen(url)
        if weburl.getcode() == 200:
            data = json.loads(weburl.read())
        for i in data["features"]:
            mag, place = i["properties"]["mag"],i["properties"]["place"]
            x,y = float(i["geometry"]["coordinates"][0]),float(i["geometry"]["coordinates"][1])
            point = Point(x,y)
            output.write({'properties':{'Place': place, 'Magnitude': mag},
            'geometry': mapping(point)})

shpname = "C:\Users\Ralph\Desktop\quakes.shp"
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_month.geojson"
GeoJSONToFC(shpname,url)

the crs variable contains the 4326 CS or WGS 1984. the error occurs on line 9 when I insert the crs variable containing the coordinate system into the collection.

Traceback (most recent call last):
  File "C:\Users\Ralph\Desktop\GeoJSONToShapefile.py", line 23, in <module>
    GeoJSONToFC(shpname,url)
  File "C:\Users\Ralph\Desktop\GeoJSONToShapefile.py", line 9, in GeoJSONToFC
    with collection(shpname, "w", crs, "ESRI Shapefile", schema) as output:
  File "C:\Python27\ArcGIS10.1\lib\site-packages\fiona\__init__.py", line 180, in open
    this_schema = schema.copy()
AttributeError: 'str' object has no attribute 'copy'

the script runs fine when I remove the crs variable but the shapefile has no projection

I have seen this post Using Fiona to write a new shapefile from scratch where the answer to defining a CS is

with fiona.open(shapeout, 'w',crs=from_epsg(3857),driver='ESRI Shapefile', schema=yourschema) as output:

I have tried this method and it gives me the exact same error as above. I also checked out the doc pages on fiona to no avail http://toblerity.org/fiona/fiona.html#module-fiona.crs


texts: OpenLayers: Feature not showing when added manually

texts:
Here my sample code for adding feature to map, but feature is not shown on map and there is no script error also. Using OpenLayers version 5.3.3.
import "ol/ol.css";
import Map from "ol/Map";
import View from "ol/View";
import WKT from "ol/format/WKT";
//import Projection from 'ol/proj/Projection';
import {fromLonLat, get as getProjection} from 'ol/proj';
import { OSM, Vector as VectorSource } from "ol/source";
import { Tile as TileLayer, Vector as VectorLayer } from "ol/layer";

var raster = new TileLayer({
  source: new OSM()
});

var wkt =
  "POLYGON ((77.264333799000042 28.705835345000025, 77.264238837000025 28.705819646000066, 77.264253835000034 28.705667240000025, 77.264259793000065 28.705570296000076, 77.264386148000085 28.705606233000026,77.264418138000053 28.705622619000053, 77.264333799000042 28.705835345000025))";

var format = new WKT();

//var feature = format.readFeature(wkt);
//feature.getGeometry().transform("EPSG:4326", "EPSG:4326");

var feature = format.readFeature(wkt, {
  dataProjection: "EPSG:4326",
  featureProjection: "EPSG:3857"
});

var geom = format.writeGeometry(feature.getGeometry());
console.info(geom);

var extent = feature.getGeometry();

var vector = new VectorLayer({
  source: new VectorSource({
    features: [feature]
  }),
  type: 'Polygon',
  projection: "EPSG:3857",
  visibility: true
});

var map = new Map({
  //layers: [raster, vector],
  target: "map",
  view: new View({
    //units: "m",
    //projection: SRS,
    //projection: "EPSG:3857",
    //center: fromLonLat([77.26433379900004, 28.705835345000025]),
    maxZoom: 12,
    zoom: 5
  })
});

map.addLayer(raster);
map.addLayer(vector);

map.getView().fit(extent, map.getSize());

//var geom1 = format.writeGeometry(map.getGeometry());
//console.info(map.projection);

Any suggestion what am I doing wrong here?


texts: ArcGIS - Geoprocessing creating temporary files (e.g. "g_g193"). What produces these and how can they be cleaned up?

texts:
I have a couple of ArcGIS extensions that are running a lot of geoprocessing models with a mix of ArcObjects. I've noticed that the some, as of yet unidentified, geoprocessing steps are creating temporary raster files (e.g. g_g193, g_g983, g_g892). The funny thing is, these files are not being directed to my selected temporary workspace or the default Arc workspace. They are actually being sent to the 'bin' folder of my extension, which is obviously an issue since these temp files should not be stored with my included DLLs and EXE files.

It's important to mention that my extension makes many calls to Arc geoprocessing, so it's difficult to location the source of these files. The interesting thing is, if I walk through all of my code using break points in .NET these files don't always get created. It's like Arc is not cleaning up these temp files if the processing is moving very quickly, but it will clean them up if given time (e.g. user steps through slowly). Regardless, I don't want these files anywhere near my 'bin' folder.

Any insight would be greatly appreciated!


texts: display multiple vector layers from single PostGis table in GeoServer using OpenLayers -- REWORDED

texts:
Prompting suggestions I’m going to attempt to reword my question

I have a PostGis table that contains overlapping vector polygons. I need to display a selection in Geoserver and provide functionality for the user to toggle on & off unwanted polygons.

I need to add multiple layers from a PostGis Table.

I’ve made a start with the code below but I’m stuck with the following:

Why are all the layers Base Layers? Land should be the only base
layer the rest overlays?
Can all overlay layers be on at once?
How do I move the layerswitcher control outside the map area? When I
click inside the map how do I show the attribution?
How do I create a new layer by script either with a polygon given a
co-ordinate string, or a point and radius?
How do I create a list of all the selected layers?

       body {
            margin: 1em;
        }
        #map {
            width: 800px;
            height: 475px;
            border: 1px solid black;
        }

        OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3;

    var map;
    function init(){
        map = new OpenLayers.Map('map');
        var Base = new OpenLayers.Layer.WMS(
            "Land",
            "http://localhost:8080/geoserver/cite/wms",
            {LAYERS: 'cite:land',
             format: 'image/png',
             isBaseLayer: true

            }
        );

        var L10299 = new OpenLayers.Layer.WMS(
            "10299",
            "http://localhost:8080/geoserver/cite/wms",
            {
            VIEWPARAMS: 'kathy:10299',
            LAYERS: 'cite:Q_Results',
            isBaseLayer: false
            }
        );

        var L43352 = new OpenLayers.Layer.WMS(
            "43352",
            "http://localhost:8080/geoserver/cite/wms",
            {
            VIEWPARAMS: 'kathy:43352',
            LAYERS: 'cite:Q_Results',
            isBaseLayer: false
            }
        );

        var L25177 = new OpenLayers.Layer.WMS(
            "25177",
            "http://localhost:8080/geoserver/cite/wms",
            {
            VIEWPARAMS: 'kathy:25177',
            LAYERS: 'cite:Q_Results',
            isBaseLayer: false
            }
        );

        map.addLayers([Base, L10299, L43352, L25177]);

        map.addControl(new OpenLayers.Control.LayerSwitcher());

        map.zoomToExtent(
            new OpenLayers.Bounds(12.996, 52.172,13.931, 52.958)
        );
    }
</script>

    OpenLayers PostGis Layers


texts: GDAL/OGR ogr.Layer.Intersection fails when producing mixed geometry results

texts:
I am trying to intersect layers using the ogr.Layer.Intersection method, in Python. This works fine for pairs of polygons, but I have found intersecting two linestring layers containing matching points can return a combination of linestrings and points, producing an invalid layer for many single-geometry type data sources e.g. Shapefile. How can I avoid the non-matching geometries without going down to the Feature level and filtering them myself?

Simplified example below, which demonstrates a POINT geometry is returned from the two linestrings.

from osgeo import ogr, osr

line1 = "LINESTRING ({0} {0}, {1} {0})".format(0, 1)
line2 = "LINESTRING ({1} {0}, {2} {0})".format(0, 1, 2)

sr = osr.SpatialReference()
sr.ImportFromEPSG(27700)

driver = ogr.GetDriverByName('MEMORY')
driver2 = ogr.GetDriverByName('ESRI Shapefile')

ds1 = driver.CreateDataSource('ds1')
ds2 = driver.CreateDataSource('ds2')
ds3 = driver.CreateDataSource('path to shapefile.shp')

layer1 = ds1.CreateLayer('layer1', sr, ogr.wkbLineString)
layer2 = ds2.CreateLayer('layer2', sr, ogr.wkbLineString)
layer3 = ds3.CreateLayer('layer3', sr, ogr.wkbLineString)

geom1 = ogr.CreateGeometryFromWkt(line1)
feat1 = ogr.Feature(layer1.GetLayerDefn())
feat1.SetGeometry(geom1)

geom2 = ogr.CreateGeometryFromWkt(line2)
feat2 = ogr.Feature(layer2.GetLayerDefn())
feat2.SetGeometry(geom2)

layer1.CreateFeature(feat1)
layer2.CreateFeature(feat2)

layer1.Intersection(layer2, layer3) # This fails, as a point is produced.

UPDATE: The exact error message received if ds3 is referring to an ESRI Shapefile on disk is "ERROR 1: Attempt to write non-linestring (POINT) geometry to ARC type shapefile.". So my main question is: how to avoid the POINT feature being produced at all, or to filter them before the layer is written to the DataSource?

UPDATE #2: Just to be clear, in my real scenario, I am not wanting to intersect two adjacent linestrings, I actually want to clip a roads Shapefile using a polygon of my study area. Coincidentally, it seems one of the road vertex points and the edge of my polygon coincide, producing the error. It worked for all the other layers I've clipped. I do not want the non-matching geometries, and believe the Intersection method should have an option to drop them if the layer has been specified as a particular type.

UPDATE #3: As a working example, consider arcpy.Intersect_analysis:

arcpy.Intersect_analysis(['shape1.shp', 'shape2.shp'], 'shape3.shp')

Even with features that intersect at a point, no multi-geometry outputs are being produced here, they must be filtered out or not produced at all - how can I get the same behavior from Python GDAL/OGR (not shapely or a database) at an ogr.Layer level?


texts: How to access all mxd files in a folder?

texts:
I'm a beginner to Python and am trying to figure out how I can access several mxds in a folder. For example, when I want to use the AnalyzeforMSD function in ArcGIS:

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Project\ReadyForMSD.mxd")
analysis = arcpy.mapping.AnalyzeForMSD(mxd)

How can I define a folder with mxds as a path, e.g. C:\Projects\data, and analyze all the mxds in it?

import os, arcpy

directory = r"C:\Project"
 for root, dirs, files in os.walk(directory):
   for myFile in files:
        fileExt = os.path.splitext(myFile)[1]
        if (fileExt == ".mxd"):
            fullPath = os.path.join(root, myFile)
            print myFile

            myMap = arcpy.mapping.MapDocument(fullPath)
            analysis = arcpy.mapping.AnalyzeForMSD(myMap)

 for key in ('messages', 'warnings', 'errors'):
 print "----" + key.upper() + "---"
 vars = analysis[key]
 for ((message, code), layerlist) in vars.iteritems():
print "    ", message, " (CODE %i)" % code
print "       applies to:",
for layer in layerlist:
    print layer.name,
print


texts: Calculating raster mean values for differnt class of another raster using ArcGIS Desktop

texts:
I have two rasters, First is LULC and another is NDVI. now  I want for each land-use class

How can we calculate avg values of NDVI in ArcGIS Desktop?

Let's say classes are Built-up, Forest, Agri, etc...

I have tried with vectorization of LULC raster and then adding the values of NDVI raster to a polygon of LULC, BUt the results are not accurate.


texts: Map data with gradient colors not showing correctly online, but works fine in ArcGIS Pro

texts: My map with gradient color symbology looks/acts correctly in ArcGIS Pro. When I successfully upload it into ArcGIS Online, the colors no longer match up with the data. I think this has something to do with the fact that the data is a one-to-many relationship for the 50 US States (multiple data points for each state over time): I think ArcGIS Online is not reading the polygon boundaries correctly because there are overlapping, existing polygons for each state(for example: Each row of data for Alabama contains the same polygon column information). Is there a way for me to display this correctly on ArcGIS Online?


texts: How to run QGIS 2.18.15 with GRASS 6.4.3

texts:
Algorithm wanted by me: v.segment is not working in QGIS 2.18.15 with grass7 https://issues.qgis.org/issues/17850, so I have to use grass6. It was working in my QGIS 2.18.1, but I had to change to 2.18.15 and I don't remember previous settings.

Starting from scratch. I install QGIS and then try to set paths to GRASS 6 from older QGIS in Processing->Options, but this is no working :(

EDIT:
Here is the option window (sorry for polish):


texts: Looking for a world database with locations of medical centers, schools, kindergartens, etc

texts: Is there some database with official information about medical centers, schools, POI's and their geo coordinates for all countries in the world?  It doesn't have to be free, but I need the GEO coordinates as well, because I need to connect them to my products!


texts:
I work in an Architectural practice concentrating on large scale urban planning and Urban Design Schemes.

Our current software is not really geared for this sort of work.

We need to be able to accurately plot building block models onto accurate cadastral information.

Is there any type of software that would allow us to do this sort of work?


texts:
I am working for a client and we chose to use OpenGeo Suite as our system. Our client will be flying a plane over farm land capturing images that we need to put in our system. From my understanding those images will be raster files. (.jpeg- just from a camera)

He will also be having GPS locations with each of the images, but in a separate file.

I can't figure out how to convert those raster files so that I can upload them into the system and be able to view them on an interactive map, in the correct spot.


texts: ArcPad Form Compatibility with previous versions?

texts:
Are forms created in ArcPad Studio 10 compatible with previous editions of ArcPad?

I have lots of units loaded with ArcPad 7 and no way to get custom forms on them (I lack ArcPad Studio or Application Builder 7). I'm trying to get basic things like pull down menus and checkboxes.


texts: How can I change the priority for MapService Process?

texts:
I've got an asynch GP service that exports mxds to pdfs.  Typically the mxd contains mapservice layers.  The same MapServices are also used by interactive clients where draw speed is important.  Drawspeed is not that important for the GP service though.

When the GP service does an export there is a noticeable spike in CPU usage on the machine hosting the MapServices.

I would like to reduce the impact on the CPU posed by this GP service which is hosted on a different machine.

To do that I'm considering creating some new MapServices specifically for use by the GP service.  I'd like for them to serve the same mxd/msd's but run at a lower priority.

Is there a way to configure this?

If not, would setting Process.PriorityClass to BelowNormal with a Server Object Extension (SOE) when IServerObjectExtension.Init gets called have any ill side effects?


texts: Choosing layers to dump in legend for a standalone application

texts:
I'm writing a standalone application using QGIS API to export maps automatically as pdf.

When I'm building my QgsComposition, I have an issue:

I have already built a list with the layers I want to show on the map (not all the layers from the original project), but my QgsComposerLegend shows all the layers.

First, I tried this:

    legend = QgsComposerLegend(c)
    legend.setTitle(u"Légende:")
    legend.model().setLayerSet(renderer.layerSet())
    [...]
    c.addComposerLegend(legend)

And the result:

But the layer set is not the one I want.

I found this post: Filter MapComposer legend to only show layers visible on the map

I adapted its code to mine and I have this:

With the code:

    legend = QgsComposerLegend(c)
    legend.setTitle(u"Légende:")
    layerGroup = QgsLayerTreeGroup()
    for layer in layers:
        layerGroup.insertLayer(0, layer)
    legend.modelV2().setRootGroup(layerGroup)
    [...]
    c.addComposerLegend(legend)

The issue is that I would keep the initial tree organization, and just suppress the layers I don't want.

I tried using the QgsLayerTreeGroup methods this way:

    legend = QgsComposerLegend(c)
    legend.setTitle(u"Légende:")
    layersToRemove = []
    for layer in legend.modelV2().rootGroup().findLayers():
        if layer.layer() not in layers:
            layersToRemove.append(layer.layer())
    for layer in layerToRemove:
        # print(layer.id())
        legend.modelV2().rootGroup().removeLayer(layer)
    [...]
    c.addComposerLegend(legend)

But I have the same result as the first code... I don't understand why, but the layers are not removed from the tree... Any idea ?


texts: Editing in ArcGIS Pro: Tracing to Modify Polygon

texts:
Currently trying to update city area polygons based on parcel lines. Some of the pre-existing area annexation polygons need to be drawn exactly coincident to the parcel line.

I normally do this in arcmap using the tracing tool, but Pro does not offer tracing while modifying a polygon, so I'll have to go back and do extra topology correction. Is there a workflow that does the equivalent?


texts: Select By Attribute: Between Two Numbers

texts:
I am trying to select house numbers that are between 200 and 499. I've been using:

"HOUSE_" >= '200' AND "HOUSE_" <= '499'

... and it's coming up with almost everything else. The numbers it does grab don't make sense at all either. I tried using a BETWEEN operator but that just gives me an error.


texts:
I am trying to mosaic very large rasters in R.  My code, taken from this helpful answer, worked perfectly for 3 test rasters (50 MB, 50 MB, 120 MB), but when I tried a larger test on 7 rasters (ranging from 1 MB to 400 MB) I returned an error.

setwd("~/raster_tiles")

gdalbuildvrt(gdalfile = "*.tif", # uses all tiffs in the current folder
             output.vrt = "test.vrt")

gdal_translate(src_dataset = "test.vrt",
               dst_dataset = "test.tif",
               output_Raster = TRUE, # returns the raster as Raster*Object
               # if TRUE, you should consider to assign
               # the whole function to an object like dem <- gddal_tr..
               options = c("BIGTIFF=YES", "COMPRESSION=LZW"))

Error in .local(.Object, ...) :

In addition: Warning message:
In system(cmd, intern = TRUE) :
  running command '"C:\Program Files\QGIS 2.18\bin\gdal_translate.exe" -of "GTiff" "geom1to7test.vrt" "geom1to7test.tif"' had status 1
Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) :
  Cannot create a RasterLayer object from this file. (file does not exist)

I can see a .vrt file in my folder of ~19 MB, but can't seem to return a GeoTIFF. Does it only have to do with the file size when I use gdal_translate on more and bigger rasters?  Any suggestions?


texts: How to make some constraint to be evaluated first in mysql?

texts:
SELECT BusinessID
FROM
  tableauxiliary
WHERE
  MBRContains(
        GeomFromText (
            'MULTIPOINT(-6.1928749092968 106.7651250907,-6.1771250907032 106.7808749093)'
            ),
            Latlong)
  AND MATCH (FullTextSearch) AGAINST ('res*' IN BOOLEAN MODE)
AND Prominent >15

I want MBRContains restriction to be executed first because it's far more restrictive. But mysql optimizer uses the FullTextSearch first.

I ended up using subquery to force MBRContains to be executed first. I want easier way.


texts: I was wondering if there is a size, or count limitation with gdaltindex. I ran the operation from within OSGeo4W via a Python script, and from within QGIS 2.2 (Windows) and had issues processing a directory with 880 GeoTiffs. The file size of the tiffs are ~ 205MB each. However, when I processed a sub-sample (four tiffs) the operation ran perfect, from within QGIS and OSGeo4W.


texts: Seeking Spatial ETL or some other tools?

texts:
I need to read geographic data from an ORACLE spatial database in Gauss Boaga SRS, convert them in WSG84 (EPSG 32632), e put them in another ORACLE spatial data base.

I've tried to use GeoKettle 2.5 but I think that there is some bugs because it seems not set SRID on the new geometries.

I need a tool that works on server side and I would not to extract data from the source database, convert them in some intermediate format, to apply the srs conversion and then put the data in the new destination database.

Is there any other open source ETL or some other tool?


texts: Creating multiple polygons within one geometry generator using QGIS?

texts: Is it possible to create multiple geom_from_wkt() elements like Using QGIS Geometry Generator to get rectangle from point? in a single geometry generator symbol layer?


texts: I have two polygon layers, one is about different land owners boundaries and is used for atlas generation, and the other is about forests inside each property. In atlas generation sometimes i want to zoom to the intersection area between the boundary layer (used for atlas iteration) with the forest layer. That is because some times a huge property can have a small part of forest so i want atlas to zoom in to that part of the property.


texts: Wait for user selection from user interface

texts:
I am making a custom maptool for QGIS
I am overriding canvasReleaseEvent.
when the user clicks on part of map canvas I call another function which identifies and creates a  list of identified feature (there is some logic to create the list).
Then this list is used to  create the selection interface like the one below:

So far everything works well. The problem is,my program doesn't wait for the result from Selection Dialog box :

    def canvasReleaseEvent(self, event):
        parrent=self.select_parent_feature('structurepoint',event)
        #here the program should wait till the result from parrent is ready
        self.main_dlg.ParrentID_txt.setText(parrent)

Here is the function which creates the selection list, shows dlg
and lets user to select

    def select_parent_feature(self,layer_name,e):
    '''This Functtion will show a list of all features in clicked point and returns user selection
    :param Layer_name: Qgis Layer name should be Already open.
    :type Layer_name: string

    :param e: is mouse event
    :type e: event
    '''

    #Start of creating list and identifieng the features
    global selection_list_header #this is the list of headers of table
    selection=""
    selection_done_flag=False
    identified_parent_list=[]
    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name() == layer_name:LayerID = lyr
    features = QgsMapToolIdentify(self.canvas).identify(e.x(), e.y(), [LayerID], QgsMapToolIdentify.TopDownAll)

    if len(features) > 0:#to make sure there is feature in clicked area
        if len(features)==1:#if there is only one item there is no need for the selection list
            feature = features[0].mFeature# the list only has one element
            selection=feature['id']
            selection_done_flag=True

        else: # show list to the user to select from
            for obj in features:
                feature = obj.mFeature
                info=[feature['id'],feature['size']]
                identified_parent_list.append(info)
                info=""

            table=self.selectionList_dlg.Selectionlist_tbl
            table.setColumnCount(len(selection_list_header))
            for i in range(len(identified_parent_list)):
                table.insertRow(i)

            column_index=0#the following line sets the headers
            for item in selection_list_header:
                table.setHorizontalHeaderItem(column_index, QTableWidgetItem(item))
                column_index+=1

            row_number=0#the following loop populates the value in table and lock the cells
            for item in identified_parent_list:
                for column_index in range(len(selection_list_header)):
                    cell_value=QTableWidgetItem(item[column_index])
                    cell_value.setFlags(QtCore.Qt.ItemIsEnabled)#lock the cell
                    table.setItem(row_number,column_index,cell_value)
                row_number+=1
            #end of logic for createn selection list DLG
            #this it where my program doesnt work
            self.selectionList_dlg.show()
            result = self.main_dlg.exec_()
           # we need to wait for the result to be ready
            if result:
                current_row_number=table.currentRow()
                if current_row_number<>-1:
                    selection=table.item(current_row_number,0).text()#we asssumed the id is located on column 0
                    selection_done_flag=True

    #I tried to  use wait until and it explodes the program
     #while not selection_done_flag:#Wait till the user select from the list
         #time.sleep(0.0001)

    return selection

As the program shows, I defined a flag to check and wait until the flag is true but that explodes the program.

the problem is the last part of the program where we show the selection list and the user has to select before the program can continue


texts: Unable to georeference Arrowsmith's "The World, on Mercators Projection" historical map

texts:
I've unsuccessfully tried to georeference this world map for several days now and am at my wits' end.

Whatever I try, however many control points I set, I get garbage when I view the resulting GeoTIFF.  I've tried https://www.georeferencer.com/ (in which I've set 60 control points with minimal deviation) and QGIS 3 (which has yielded nothing but "space-warped art") and I have simply no luck.
Is there an issue with the map, maybe, anyone could identify?  Is there any chance some kind person could have a look what might be the issue?
Map overview: https://www.davidrumsey.com/luna/servlet/detail/RUMSEY~8~1~2722~260002:The-World,-on-Mercators-Projection-
Download link: https://www.davidrumsey.com/rumsey/download.pl?image=/D0026/4613002.sid

Bottom control point is Cape Town, is that just a visual error or a real one?

Finished product, there's actually a thin blue line where the base map shines through in the middle that I don't know how to get rid of, but other than that it looks great.  Thanks again!


texts: Can any one tell me how to convert line features into polygon features in ArcGIS.  I am using ArcGIS 10. I searched the arc toolbox for any tools for this conversion, but I couldn't find any.  Please help me regarding this.


texts:
I am trying to make a 1km square grid using QGIS 2.18.26.

I am using EPSG:4326. I need it to be exactly 1km square because that's how the variable data has been collected.

I am struggling to get the grid to be exactly 1km square, I have attempted using Vector Create Grid, and MMQGIS create grid layer, but I am unsure how to set the grid units to 1km square.


texts: Update ArcPy version to avoid AttributeError: 'module' object has no attribute 'GetClassNames'?

texts:
This is probably the answer to my previous question: Listing Elements or Datasets of Layer gives ERROR 000840 value is not Feature Layer?.

The version of arcpy that I am using does not include GetNAClassNames

  naClasses = arcpy.na.GetNAClassNames(serviceAreaLayer)

AttributeError: 'module' object has no attribute 'GetClassNames'

According to this help topic, GetNAClassNames is in ArcGIS 10.1. Why is it not in my version?


texts:
I have a problem getting a large number of facilities and incidents to draw upon solving.

I have about 380 facilities and incidents, that I am attempting to find routes between.

I need the actual network geometry, so and OD Cost Matrix won't work. What happens is I load in everything and attempt to solve and it will get to the very end (99%), and then say it is completed, but there are no actual routes. I have got it to route to 299 incidents and 299 facilities, but with 308, the solve fails.

I began creating a model to run it outside of ArcMap, because I thought that too many routes trying to draw might be the issue, but the model is returning an error.

I have attached a photo of the model.

Any suggestions on how to fix this? Or what the problem might be? Is the memory just insufficient or something?

I'm using ArcGIS 10.0, with the Network Analyst Extension.


texts: Transparent color of Selection in Attribute Table

texts: I want to check the numbers of the selected vector in the attribute table. In my setting, the following is my table after I selected some vectors. As you can see, I am not able to look at any values. What setting should I need to resolve this issue?


texts: Sudden errors intersecting NDVI raster data with overlapping polygons in GME

texts:
I am experiencing errors with Geospatial Modeling Environment (GME) that suddenly developed after things were working perfectly.

Description of software installed: ArcGIS 10.3.1 for Desktop and ArcGIS Service Pack: 0 (build 0); R version 3.3.2; Geospatial Modelling Environment (Version 0.7.4.0)

What I'm trying to do

I am using the tool "isectpolyrst." My polygon file has 5718 rows
with overlapping buffers. Buffer size is 90m around points. I am
trying to extract Normalized Difference Vegetation Index raster data
(NDVI; from MODIS) to each polygon. NDVI values range from -2000 to 9000
I do not get an error message. However, the values of the raster I am trying to extract do not match the output. See below.
The tool processes very, very quickly. The output value for each
polygon is identical (mean = -1.797693e+308). These values are not
correct (ie there is variation and the values do not match) The
tricky part is the tool worked properly until a few weeks ago. I did
not upgrade GIS, nor did I change the files at all.

3. What I've tried to trouble shoot:

Double check projections are the same for polygon and raster data.
Downloading NDVI files again. Perhaps they got corrupted?
Scaling my data to 0-100 values (subtract by minimum value, divide by range, multiply by 100).
Converting to different raster file type (TIFF, IMAGINE, GRID)
Assigning attribute table to NDVI raster data (no luck)
Try on different machines (all have Windows)
Using the Zonal Stats as Table II in ArcMap. This is a finicky tool and I can never get it to work on >300 rows.

I AM STUMPED!!! Does anyone have any ideas on how to proceed? This is an issue I've heard from other folks, but no one found a work around in GME.


texts: WCS 1.1.1 OGC Complaint Test Failed in Team Engine

texts:
I have been testing my GeoServer 2.12.1 for OGC Compliant to WCS 1.1.1 Standard. But I got following errors while testing with the OGC Team Engine.

How do I resolve the issues?

wcs:GetCoverage_RangeSubset_InterpolationMethod_Missing - Failed

  Request d1e19427_1:    Method: POST    URL:
  http://BASEURL/geoserver/wcs?    Body:

<GetCoverage xmlns="http://www.opengis.net/wcs/1.1.1" xmlns:gml="http://www.opengis.net/gml"
             xmlns:ows11="http://www.opengis.net/ows/1.1"
             xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
             xmlns:ctl="http://www.occamlab.com/ctl"
             xmlns:ctlp="http://www.occamlab.com/te/parsers"
             xmlns:fn="http://www.w3.org/2005/02/xpath-functions"
             xmlns:gen="java:com.occamlab.te.Generator"
             xmlns:html="http://www.w3.org/1999/xhtml"
             xmlns:ows="http://www.opengis.net/ows/1.1"
             xmlns:saxon="http://saxon.sf.net/"
             xmlns:te="http://www.occamlab.com/te"
             xmlns:tec="java:com.occamlab.te.TECore"
             xmlns:tems="java:com.occamlab.te.web.MonitorServlet"
             xmlns:tep="http://teamengine.sourceforge.net/parsers"
             xmlns:wcs="http://schemas.opengis.net/wcs/1.1.1"
             xmlns:xi="http://www.w3.org/2001/XInclude"
             xmlns:xlink="http://www.w3.org/1999/xlink"
             xmlns:xs="http://www.w3.org/2001/XMLSchema"
             xmlns:xsd="http://www.w3.org/2001/XMLSchema"
             service="WCS"
             version="1.1.1"
             xsi:schemaLocation="http://www.opengis.net/wcs/1.1.1 http://schemas.opengis.net/wcs/1.1.1/wcsAll.xsd">
   <ows11:Identifier>topp:rotated</ows11:Identifier>
   <DomainSubset>
      <ows11:BoundingBox crs="urn:ogc:def:crs:OGC:1.3:CRS84">
         <ows11:LowerCorner>7.7634301664746515 45.14713380418506</ows11:LowerCorner>
         <ows11:UpperCorner>7.764350661575157 45.14763319238466</ows11:UpperCorner>
      </ows11:BoundingBox>
   </DomainSubset>
   <RangeSubset>
      <FieldSubset>
         <ows11:Identifier>contents</ows11:Identifier>
         <InterpolationType/>
      </FieldSubset>
   </RangeSubset>
   <Output format="image/tiff"/>
</GetCoverage>

  Error in call to extension function {public org.w3c.dom.NodeList
  com.occamlab.te.TECore.request(org.w3c.dom.Document,java.lang.String)
  throws java.lang.Throwable}: Exception in extension function
  java.lang.reflect.InvocationTargetException Result: Failed

wcs:GetCoverage_Store_Bogus - Failed

  Request d1e24800_1:    Method: GET    URL:
  http://BASEURL/geoserver/wcs?&service=WCS&request=GetCoverage&version=1.1.1&identifier=topp:rotated&BoundingBox=7.7634301664746515,45.14713380418506,7.764350661575157,45.14763319238466,urn:ogc:def:crs:OGC:1.3:CRS84&Format=image/tiff&Store=bogus
  Error in call to extension function {public org.w3c.dom.NodeList
  com.occamlab.te.TECore.request(org.w3c.dom.Document,java.lang.String)
  throws java.lang.Throwable}: Exception in extension function
  java.lang.reflect.InvocationTargetException
  Result: Failed

I have changed the input for bogus to various values but still not working. The procedure followed for testing the WCS is from http://docs.geoserver.org/2.12.1/developer/cite-test-guide/index.html#run-wcs-1-1-tests.

Also, the WCS Capabilties base URL to be entered for testing the WCS Compliance as mentioned in the post https://github.com/opengeospatial/ets-wps10/issues/9#issuecomment-253070197.


texts: Plotting polygons as separate plots using Python

texts:
I have a shapefile with several polygons in a geodataframe. I want to plot each individual polygon in their own plot.
I tried
for pol in z3:
  plt.figure()
  plt.plot(pol)
  plt.show()

But this just loops through the columns of the data frame.


texts: Leaflet - global variable in onEachFeature returns value of next feature

texts:
Thanks to TomazicM I got my GeoJSON import up and running. Here's the whole code, I highlight the problematic bits thereafter

    var map = L.map('map', {
        crs: L.CRS.Simple,
        minZoom: -1
    });

//TODO - set the image width and height to w and h
    var w = 1700;
    var h = 2200;

    //the stuff to be able to work with xy instead of latlon
    var yx = L.latLng;
    var xy = function(x, y) {
        if (L.Util.isArray(x)) {    // When doing xy([x, y]);
            return yx(x[1], x[0]);
        }
        return yx(y, x);  // When doing xy(x, y);
    };

    var bounds = [xy(0, 0), xy(w, h)];
        map.fitBounds(bounds).setMaxBounds(bounds); //this makes the map fit the bounds we just defined

//TODO - insert image name here
    var image = L.imageOverlay('img/kisandra.png', bounds).addTo(map);

    //create map icon class by extending L.Icon
    var TaladasIcon = L.Icon.extend({
    options: {
        iconSize:     [50, 50],
        iconAnchor:   [25, 25],
        popupAnchor:  [0, -25]
        }
    });

    //define custom icons
    var metropolisIcon = new TaladasIcon({iconUrl: '/img/metropolis.png'}),
        cityIcon = new TaladasIcon({iconUrl: 'img/city.png'}),
        townIcon = new TaladasIcon({iconUrl: 'img/town.png'}),
        villageIcon = new TaladasIcon({iconUrl: 'img/village.png'}),
        innIcon = new TaladasIcon({iconUrl: 'img/tankard.png'});

    var featureIcon     = villageIcon;

function forEachFeature(feature, layer) {
        var popupContent = "
 " + feature.properties.name +
            "
"+ feature.properties.description +
            "
"+ feature.properties.icon +"
";
        featureIcon     = window[feature.properties.icon];
        layer.bindPopup(popupContent);
    }

var jsonMetropolis = L.geoJSON([], {
  onEachFeature: forEachFeature,
  pointToLayer: function (feature, yx) {
    return L.marker(yx, {
      icon: featureIcon
    });
  }
});

$.ajax({
  dataType: "json",
  url: "data/metropolises.geojson",
  success: function(data) {
    $(data.features).each(function(key, data) {
      jsonMetropolis.addData(data);
    });
    jsonMetropolis.addTo(map);
  }
});

As you can see after I defined my custom icons, I define a variable called featureIcon and assign it a random icon, just so it's not null var featureIcon     = villageIcon;
Then in the forEachFeature function I set the variable equal to a field in the geoJSON         featureIcon     = window[feature.properties.icon]; and later the marker is told to use said icon L.marker(yx, {icon: featureIcon})

Funny thing is, while the popups that are set in forEachFeature too, are displaying correctly the property field's contents, the icon is offset by one feature. The first feature displays the pre-set icon, the second feature displays the icon specified in feature 1, the third the icon from feature 2, etc.

What am I doing wrong this time? Wrong use of the global variable thing window[] ? Or `something else?


texts: Export multiple attachments per point feature from fgdb using python script

texts:
I use the following script for export of attachments from an arcgis fgdb. It is used inside a script tool. I have multiple attachments (relID) for every pointfeature (GLobalID) and I would like to modify it to include export all of those.
How do I achieve this? I would like to be able to distinguish between  attachments which relate to the same pointfeature somehow in the output folder.
Edit:
New code:
from arcpy import da
import os

attachTable = arcpy.GetParameterAsText(0)
origTable = arcpy.GetParameterAsText(1)
nameField = arcpy.GetParameterAsText(2)
fileLocation = arcpy.GetParameterAsText(3)

originCursor = da.SearchCursor(origTable, ['GlobalID', nameField])
for origRow in originCursor:
    origID = origRow[0]
    origName = origRow[1]

    attachCursor = da.SearchCursor(attachTable, ['DATA', 'ATT_NAME', 'REL_GLOBALID'])
    for attRow in attachCursor:
        binaryRep = attRow[0]
        fileName = attRow[1]
        relID = attRow[2]
if relID == origID:
    open(fileLocation + os.sep + origName + ".jpg", 'wb').write(binaryRep.tobytes())

del originCursor


texts:
I am new to ArcGIS Desktop 10.  I am using a students version on my home laptop.  I have lost my table of contents completely from the left side of my screen.

Does anyone know where I can find it?


texts: How do I install a plugin when I have the .zip file for it

texts:
I have the zipped file toLabeledLayer-1.0.zip and want to install and use it.  I am using QGIS 2.18.6.  I have unzipped the file and copied it into .qgis2\python\plugins as suggested here but it does not appear in the plugin manager.  I have also copied to C:\Program Files\QGIS 2.18\apps\qgis\python\plugins and
I have also tried to add the source as a repository under settings in the Plugin Manager as suggested in Section 1.3 of this page.  Can anyone explain how to install this plugin to me?  Does it not work on this version of QGIS?


texts:
I am trying to convert regular python script to toolbox. This is the regular code which works:

import arcpy
strOutpath = "C:\Users\name\Desktop\New Folder"
mxd = arcpy.mapping.MapDocument("CURRENT")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
 mxd.dataDrivenPages.currentPageID = pageNum
 pageorder = mxd.dataDrivenPages.pageRow.BR
 pageindex = mxd.dataDrivenPages.pageRow.KARTA
 pdfname = pageorder + "-" + pageindex + ".pdf"
 arcpy.mapping.ExportToPDF(mxd, strOutpath + pdfname, resolution="300dpi", image_compression="Adaptive", jpeg_compression_quality="85", colorspace="RGB", picture_symbol="VECTORIZE_BITMAP")
del mxd

And this is the my example for script to input as script for toolbox:

import arcpy
strOutpath = arcpy.GetParameterAsText(0)
mxd = arcpy.GetParameterAsText(1)
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
 mxd.dataDrivenPages.currentPageID = pageNum
 pageorder = arcpy.GetParameterAsText(2)
 pageindex = arcpy.GetParameterAsText(3)
 pdfname = pageorder + "-" + pageindex + ".pdf"
 arcpy.mapping.ExportToPDF(mxd, strOutpath + pdfname, resolution="300dpi", image_compression="Adaptive", jpeg_compression_quality="85", colorspace="RGB", picture_symbol="VECTORIZE_BITMAP")

As parameters, I entered this:

strOutpath = FOLDER
mxd = ArcMap document
pageorder = field
pageindex = field

And when I try to input the values, I get empty fields. There is clearly something wrong with the second code and input variables, but I don't know what.

EDIT:

When I edit my code, based on the comment ArtWork21 gave, when only output folder and mxd are user-generated, this is the code:

import arcpy
strOutpath = arcpy.GetParameterAsText(0)
mxd = arcpy.GetParameterAsText(1)
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
 mxd.dataDrivenPages.currentPageID = pageNum
 pageorder = mxd.dataDrivenPages.pageRow.BR
 pageindex = mxd.dataDrivenPages.pageRow.KARTA
 pdfname = pageorder + "-" + pageindex + ".pdf"
 arcpy.mapping.ExportToPDF(mxd, strOutpath + '\\' + pdfname, resolution="300dpi", image_compression="Adaptive", jpeg_compression_quality="85", colorspace="RGB", picture_symbol="VECTORIZE_BITMAP")
del mxd

This is the error:

Traceback (most recent call last):
File "C:\Users\name\Desktop\script.py", line 4, in <module>
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
AttributeError: 'unicode' object has no attribute 'dataDrivenPages'
Failed to execute (Script).

And these are my variables:

So, now I now I am just trying zo get the last one script to work, where I only enter output folder, and input mxd.


texts: How to get correct area calculation using Google Mercator?

texts:
i use QGIS with openlayers-plugin and i´ve loaded google satellite as a layer. what i want to do is making a shape-layer (polygon) with the goal to get the area values of the polygons. i´ve activated on-the-fly-transformation and chosen Google Mercator (EPSG:900913) for the project and the shape-layer.

If i compare the values of the measure line tool in qgis with the one in google earth (desktop version) to make a rough estimate, i get two different values. they are same if i activate ellipsoidal on the measure line tool in QGIS. The area calculation takes places with ellipsoidal deactivated, so the area result in the attribute table is wrong. I know that the values of google earth are correct.

Does someone know a way to get an correct area calculation of the polygons?


texts:
I am new to GIS. I am done making a map, using ArcGIS Online webmap. I have successfully installed QGIS 2.14, now I want to retrieve the map and layers from ArcGIS to be used in QGIS.

What steps should I take to download the map (data) from ArcGIS online?


texts:
My goal is to reuse the default edits tool of QGIS in a custom App.
But I would like to modify some behavior with those tools.

For instance, the MapTool beside the Edit Node Action.

iface.actionNodeTool()

When I activate the tool with the GUI I can get the tool beside the QAction with:

iface.mapCanvas().mapTool()

Give me this:

<qgis._gui.QgsMapTool object at 0x0000000011167A60>

Does that mean this component inherits directly from the QgsMapTool ?
And we cannot reuse or override the logic ?

If I want the same behavior (and more) I have to rewrite all the logic ?

Last example for illustrate, with the polygon edit tool.

For terminate a feature, you have to right click. I like all the tool behavior, but I would like to terminate with a double click. Do I have to rewrite a tool for QgsMapTool or can I reuse the existing tool ?

If you have any ideas


texts: How to register multiple events on multiple layers with Openlayers

texts:
Is there a way to register multiple events on multiple objects with one command?

In the API doc one can find several solutions(http://dev.openlayers.org/docs/files/OpenLayers/Events-js.html#OpenLayers.Events.register), but for me it seems, that the combination of multiple layer and multiple events is not possible.

Conretely I would like to register the loadstart and loadend events on all of my layers which are about 15.

I tried among others:

events.on({
    "loadstart": function1,
    "loadend": function2,
    scope: [layer1, layer2, layer3]
    });

But it didn't work (it seemed that one cannot start with events.on, the object has to be declared in the front)....


texts:
I am trying to find shorthest path between multiple points. I am using new route under Network Analyst in ArcMap.

But, if there is no road between points due to lake or sea, how can I find shortest path between multiple points? For instance, in the picture below, between point(2) and point (11), there is a lake.

Perhaps, I can create line inside the lakes by using ArcMap. Do you have any suggestion how can I do this ?


texts: Discrepancy between GTiff and MEM drivers

texts:
I was expecting that both the drivers would behave the same way, except that 'MEM' will not create the file.

However, in this test 'MEM' pass, but 'GTiff' does not.
Can someone explain why?

class TestMEMVsGTiff(unittest.TestCase):

    @staticmethod
    def check(driver_type):

        temp_tif = tempfile.mktemp(suffix='.tif')
        data = np.array([[[4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 7, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 2, 2, 7.],
                         [4, 7, 7, 10, 2, 7.]],
                        [[2, 0, 0, 0, 0, 0.],
                         [2, 0, 0, 0, 0, 2.],
                         [0, 1., 0, 0, 0, 1.],
                         [0, 0, 0, 0, 0, 2],
                         [0, 0, 0, 0, 0, 0.],
                         [0, 0, 0, 0, 0, 0.]]], dtype=np.float32)
        src_ds = gdal.GetDriverByName(driver_type).Create(temp_tif, 6, 6, 2,
                                                    gdalconst.GDT_Float32)

        src_ds.GetRasterBand(1).WriteArray(data[0, :, :])
        src_ds.GetRasterBand(1).SetNoDataValue(2)
        src_ds.GetRasterBand(2).WriteArray(data[1, :, :])
        # src_ds.GetRasterBand(2).SetNoDataValue()
        src_ds.SetGeoTransform([10, 1, 0, 10, 0, -1])
        src_ds.FlushCache()

        dst_ds = gdal.GetDriverByName('MEM').Create('', 3, 3, 2,
                                                    gdalconst.GDT_Float32)
        dst_ds.GetRasterBand(1).SetNoDataValue(3)
        dst_ds.GetRasterBand(1).Fill(3)
        dst_ds.SetGeoTransform([10, 2, 0, 10, 0, -2])

        gdal.ReprojectImage(src_ds, dst_ds, '', '', gdal.GRA_Average)
        got_data = dst_ds.GetRasterBand(1).ReadAsArray()
        expected_data = np.array([[5.5, 7, 7],
                                  [5.5, 7, 7],
                                  [5.5, 8, 7]])
        np.testing.assert_array_equal(got_data, expected_data)
        band2 = dst_ds.GetRasterBand(2).ReadAsArray()
        np.testing.assert_array_equal(band2, np.array([[1., 0., 0.5],
                                                       [0.25, 0., 0.75],
                                                       [0., 0., 0.]]))

    def test_mem(self):
        self.check('MEM')

    def test_gtiff(self):
        self.check('GTiff')


texts: Memory Layer features not displayed after deleting some features

texts:
I am developing a plugin for QGIS and have hit a snag that does not seem to have happened before.

After having created a point memory layer with attributes, I need to delete features. The following code was used to delete same, however, after this operation, all features disappear from the mapcanvas. The features are displayed in the browser, with the exception of the deleted ones, and the layer is still shown as an existing layer in the layer legend.

Here's snippets of the code:

def createTmpLyr(self):
    uri = "Point?crs=epsg:4248&field=id:integer&field=Pt:string(20)&index=yes"
    self.vl = QgsVectorLayer(uri, "temporary_points", "memory")
    QgsMapLayerRegistry.instance().addMapLayer(self.vl)

more code that adds point features successfully...

Code that generates a selection of points to be deleted

    self.vl.setSelectedFeatures(listOfPts)
    ntBndPts=self.vl.selectedFeatures()

    self.delNtBndPts(ntBndPts)
    self.vl.commitChanges()

def delNtBndPts(self, pts):
    print pts
    with edit(self.vl):
        for feature in pts:
            self.vl.deleteFeatures([feature.id()])


texts: MapServer: apply a symbol to certain data?

texts:
I have a shp file which contains some points. Each point has a "POINTTYPE" attribute which contains text.

Is it possible to apply a different symbol basing on this attribute? For example, if pointtype is "streetlamp" use a square icon, if it is "tree" use another one, etc...?

Thanks

EDIT:

Here is the final mapfile

LAYER
    NAME bla
    TYPE POINT
    GROUP bla
    DATA "bla.shp"
    STATUS OFF
    CLASS
      EXPRESSION ('[pointtype]' = 'streetlamp')
      STYLE
        SYMBOL 'circle'
        SIZE 8
        COLOR 255 0 170
        OUTLINECOLOR 0 0 0
        ANTIALIAS TRUE
      END
    END
    # Other classes can go here, they will cumulate
END


texts:
When I try to download Aster level 1T data from  the USGS Earth Explorer site it asks me to enter my username and password over and over when I click on download link of my desired product.

What should I do?


texts: Passing WorkspaceName to another thread for opening?

texts:
My winforms app gets an IWorkspaceName from an IDatasetName via the GxObject selected in a GxDialog. I want to pass it across to a worker thread for opening.

I know how to write multithreaded arcobjects code - particularly the need to create the workspace factory class on the worker thread but my app makes no assumptions as to the type of workspace factory. I've considered passing it as a serialized xml string, but the documentation seems strangely vague on this topic. Does anyone have code to do this?

On preview,

I'm asking much the same as this question. Has anyone cracked the xml serialization / deserialization? It seems the simplest way to go or is there a better way?


texts:
I'm trying to learn about the SAR image formation process.  From reference such as Carrara's Spotlight Synthetic Aperture Radar book, it appears that after various compensations, a direct 2D fft can be applied to 'compress' the data.

I understand range is in 'slow time', and crossrange/azimuth is in 'fast time'.
 However, wouldn't applying the 2D fft directly result in the frequency domain?
 What are the units in range/crossrange for raw, unfocused SAR data that allow us to directly apply a 2D fft to get a 'compressed' image (I believe this would result in a single look complex image?).


texts: How to get values of each cell in raster attribute table?

texts: How do I get the values of each cell in the raster attribute table. Usually the attribute table gives you the value and number of counts. But, in my case I want each cell listed and corresponding cell values.


texts: Convert Lat/Lon sexagesimal degrees to UTM coordinates

texts:
I'm using proj4.js to convert geographic coordinates (degrees, minutes and seconds) to UTM coordinates (ETRS89 Zone 30). When I get my coordinates result of the transformation I'm checking them with several web coordinates converters and I never get the same result (I mean really different coordinates) so I don't know if I'm doing something wrong or what.

To make the conversion, I first transform the coordinates from degrees, minutes, seconds to decimal degrees and then make the conversion like this:

var firstProjection = "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs";
this.secondProjection = "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs";
var longLat = this.toDecimalsDegrees(); //Array with coordinates in decimal degrees

var result = proj4(firstProjection, this.secondProjection, [longLat[0], longLat[1]]);

Finally, it is correct the longitude coordinate first and latitude then?

The coordinates I'm using are:

Longitude: 6º 2' 32.929''
Latitude: 37º 21' 41.131''

And once I have converted them to decimal degrees:

Longitude: 6.042480466276646
Latitude: 37.361425499749146

The result of the reprojection that I get with the code mentiones above is:

X: 1301616.242
Y: 4173536.779

Meanwhile in other web conversors the rresult is:

X: 761940.319
Y: 4139140.14


texts: Extracting a circle from a raster using GDAL

texts:
I can extract a rectangle from a raster using gdalwarp:

gdalwarp -te <x_min> <y_min> <x_max> <y_max> input.tif clipped_output.tif

Is there a way to do this in GDAL using a circle instead of a rectangle?


texts: Open a mutliband raster, edit values in one of the bands, and overwrite the raster using Python/Rasterio?

texts:
I have numerous multiband rasters (36 band) in which preserving the band order is important for future processing. One of the bands often contains erroneous data and I'd like to be able to change all the values in within that band to 0, but keep the band in place. Then ideally overwrite the changes to the existing raster. Is there a way to do this without having to open the band in question, make the change to 0, and then re-stack all the other bands along with the edited band into a new file? Or at the very least save the changes to a new raster without having to call up the other 35 bands individually into the stacking function?

For instance

import rasterio as rio
import numpy as np

inraster = "path_to_a_mutliband_raster"
rast = rio.open(inraster)
badband = rast.read(7)
badbad[badband != 0] = 0

From there I'm not sure how to save those changes back into the mutliband raster.


texts: Releasing user-schema locks when updating user-schema

texts:
Putting together a scripted event where we update user-schema (prod) in one sde instance from a dev user-schema instance.  Both sdes are Oracle 12c.  Replication is not used and the destination objects (6 feature classes and 15 tables) are not versioned.  We want to release all locks so that we can delete the objects from our production SDE and copy the same objects from the dev SDE.

Using the admin connection, can I release locks/services at the user-schema level or does it have to be performed at the master level? If so, is the cadence of events below correct?

import os
import arcpy
import datetime

start = datetime.datetime.now()

print 'start run: %s\n' % (start)

#Parameters
prodadmin = r"C:/Users/101Connections/Prod/Admin@Prod.sde" #Prod Admin
faprod = r"C:/Users/101Connections/Prod/Facilities@Prod.sde" #Prod(o)
fadev = r"C:/Users/101Connections/Dev/facilities@Dev.sde" #Dev(i)

fatabs = ["1","2","3"] #FACILITIES tables
fafcls = ["A","B","C"] #FACILITIESs feature classes

arcpy.env.workspace = prodadmin

#Release connections
users = arcpy.ListUsers()
for user in users:
    print("Username: {0}, Connected at: {1}".format(
        user.Name, user.ConnectionTime))
arcpy.DisconnectUser(faprod, "ALL")
arcpy.AcceptConnections(faprod, False)

#Change workspace for Copy+configuration
arcpy.env.workspace = faprod

#delete tables in faprod. clears out all reclass objects.
try:
    for tb in fatabs:
        if arcpy.Exists(tb):
            arcpy.Delete_management(tb)
    print "fatabs success!"

except:
    print arcpy.GetMessages(2)

#copy tables from fadev to faprod. exists to catch duplicate objects from
reclass.
try:
    for tb in fatabs:
        if arcpy.Exists(tb):
            print tb + "Already here!"
        else:
            arcpy.Copy_management(fadev + "/" + tb, faprod + "/" + tb)
            print tb + "Created!"
except:
    print arcpy.GetMessages(2)

print 'finished run: %s\n\n' % (datetime.datetime.now() - start)


texts: How to determine if a raster has rowattributes?

texts:
I have a raster file which is the result of extract by mask geoprocessing tool. If it has no rows then I take it as no overlap between input raster and shapefile. Following is the way I am reading it to determine whether a row exists or not. But it gives me zero answer even if a row exists. Please help me to better control this process.

    a = 0
    grid1 = arcpy.sa.ExtractByMask("iras",fc)
    grid1.save(optfolder + '/oras')
    arcpy.BuildRasterAttributeTable_management(optfolder + '/oras',"Overwrite")
    grid2 = arcpy.SearchCursor(optfolder + '/oras',"","","Count","")
    rows = grid2.next()
    if rows == None:
        count = 0
        a = a + count
    else:
        count = rows.getValue("Count")
        a = a + count

    del rows


texts:
I'm having a scaling problem when using DEMto3D and extracting STL's for later manipulation.

I have a merged TIFF file for Australia and am attempting to extract models of each state so I can print out. I want to print each state individually and mount them on a frame later so the scale of each model is important.

CRS for everything is GDA94.

My workflow is as follows.

Use clipper to extract state raster (using state boundary
shapofile) and input as a separate layer.
Use Demto3d on each layer to extract the states to a STL.

Whilst I can get the models saved as STL without problems, when I try to assemble them in a 3d modelling program, they are about 5-7% different in size and do not line up.

Checking each TIFF file and overlaying in GIMP shows them aligning perfectly.
DEMto3D settings are the same for each model as are the clipper settings.

There must be a projection or scaling issue but I cannot figure out, where exactly I went wrong. Nor understand why because everything I can find is set for GDA94.


texts:
I am exporting aerials with world files from data view.  I have 21 aerials I exported and each one has the credits for the service layer.  I know how to remove the credits in the layout view, but how do you do so in the data view?

I'm working in 10.2.


texts: How to select a subset of a PostGIS table when using the QGIS Offline Editing plugin?

texts: When using the Offline Editing plugin to 'Convert to an offline project' is it possible to make a selection of the dataset rather than 'offline' the whole table?


texts:
I am trying to call the setSuppress function

I want to set the suppress to 2=SuppressOff by writing in Python Console QgsEditFormConfig.setSuppress(2) but I am receiving the following error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: QgsEditFormConfig.setSuppress(QgsEditFormConfig.FeatureFormSuppress): first argument of unbound method must have type 'QgsEditFormConfig'

Does anyone know the right syntax for that the setSuppress function?


texts: SQL Condition Statement in Select By Attributes query

texts:
I am attempting to execute a Select query in Select By Attributes where the return record set I'm looking for includes all records where the String value is greater than 7 characters long.

For ex., I have records in my table where the value in the "LeftCityCo" field is 2395187, and other records where the value in the "LeftCityCo" field is 2395934ANO. I'm seeking to obtain a Selection that includes all records with values in this field greater than 7 characters long.

I've tried using the Left function and attempted the following query but returned an error:

Select * From Streets Where
"LeftCityCo" = Left( "LeftCityCo", 7 )


texts: Clip linesegment at a polygon boundary using turfjs

texts:
Having a linestring and a polygon, Is there a way to clip the linesegment at the boundary of the polygon and retain the shape of the linestring which runs inside the polygon?

Example Linestring

  {
    "type": "FeatureCollection",
    "features": [
      {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "LineString",
          "coordinates": [
            [
              -43.154296875,
              38.34165619279595
            ],
            [
              -11.42578125,
              46.13417004624326
            ]
          ]
        }
      }
    ]
  }

Example Polygon

  {
    "type": "FeatureCollection",
    "features": [
      {
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Polygon",
          "coordinates": [
            [
              [
                -35.947265625,
                28.998531814051795
              ],
              [
                -21.62109375,
                28.998531814051795
              ],
              [
                -21.62109375,
                48.3416461723746
              ],
              [
                -35.947265625,
                48.3416461723746
              ],
              [
                -35.947265625,
                28.998531814051795
              ]
            ]
          ]
        }
      }
    ]
  }

I want to know how I can clip the linesegment which crosses the polygon as a separate feature? I tried turf's lineSplit and lineIntersect. But was not able to combine multiple functions to get the required output.
Trying to understand turfs functions. Any help will be very helpful.

Thnak you.


texts:
I am having trouble getting Grass 7.2 to work in QGis 3.2. As instructed per kyngchaos, I used the following prepend PATH to get rid of error "env: python3: No such file or directory"

/Library/Frameworks/Python.framework/Versions/3.6/bin:/Library/Frameworks/GDAL.framework/Versions/2.2/Programs:

Now I am getting the following error when trying to run any Grass algorithm, with all algorithms aborting:

"WARNING: Default locale settings are missing. GRASS running with C locale."

I installed on a fresh MacOS High Sierra 13.6. Before QGIS installation I installed python 3.6 as well as GDAL as per instructions. It seems I am not the only one having this problem, but I have not been able to find any workaround.

QGIS installation information:

QGIS version
3.2.0-Bonn
QGIS code branch
Release 3.2
Compiled against Qt
5.9.3
Running against Qt
5.9.3
Compiled against GDAL/OGR
2.2.4
Running against GDAL/OGR
2.2.4
Compiled against GEOS
3.6.2-CAPI-1.10.2
Running against GEOS
3.6.2-CAPI-1.10.2 4d2925d6
PostgreSQL Client Version
9.6.6
SpatiaLite Version
4.4.0-RC0
QWT Version
6.1.3
QScintilla2 Version
2.10.3
PROJ.4 Version
493


texts: How to filter date values using OGR 2.4.0?

texts:
I am trying to replicate a working Oracle SQL query into ogr2ogr, but having issues with a date field where clause.

My query is:

select BLOCK_STATUS_DATE, GEOMETRY
from WHSE_FOREST_TENURE.FTEN_CUT_BLOCK_POLY_SVW
where BLOCK_STATUS_DATE > '01-JUL-18'
and rownum < 9
order by BLOCK_STATUS_DATE asc

This query returns a logical result, with dates that are all newer than July 1, 2018.

05-JUL-19   [MDSYS.SDO_GEOMETRY]
05-JUL-19   [MDSYS.SDO_GEOMETRY]
05-JUL-19   [MDSYS.SDO_GEOMETRY]
18-JUL-19   [MDSYS.SDO_GEOMETRY]
18-NOV-19   [MDSYS.SDO_GEOMETRY]
18-NOV-19   [MDSYS.SDO_GEOMETRY]
18-NOV-19   [MDSYS.SDO_GEOMETRY]

The same query, run through an ogr2ogr command, returns 8 rows with dates ranging from 2002 to 2019. It seems the SQL dialect is getting confused between Oracle SQL and whatever dialect ogr2ogr uses.

For completeness, here is my total ogr2ogr command:

ogr2ogr.exe -a_srs epsg:3005 -f "GPKG" T:\test\RKB_cutblocks.gpkg OCI:USER/PASSWORD@IDWPROD1:no_Table -progress -sql @T:\tempQueryFolder\query.sql -lco IDENTIFIER='layer1' -lco DESCRIPTION='layer1'

If anyone knows how I can cast the date in my query, or specify a proper SQL dialect to ogr2ogr, or any other solution, please let me know!


texts: Converting Sentinel 2 jp2 bands to TIFF with gdal

texts:
I want to convert a .jp2 file from the Sentinel 2 project to .tif with gdal. For some bands it works (e.g. B01) and for some bands it always fails (e.g. B04). I use the sentinel public image browser data (e.g this product).

When I run gdal_translate B04.jp2 B04.tif I get the following output:

  Input file size is 10980, 10980
  0maximum number of samples exceeded (120560400 > 67108864)
  error: cannot decode code stream
  ERROR 1: B04.jp2, band 1: IReadBlock failed at X offset 0, Y offset 0

My gdal version is 2.3.0, released 2018/05/04


texts:
I have two shapefiles:

a shapefile of plant communities: each polygon there covers a wide swath of land and names the plant community there.
I also have another shapefile that identifies a bunch of tree crowns and the specie of each.

Given the polygons of the second shapefile are much smaller than the first, how can I know the plant community that each crown belongs to?

I know pretty much most object oriented programming languages, R and matlab and no knowledge of ArcGIS. ENVI.


texts: Enabling Drag and Drop from TOC to Python Script Tool?

texts:
I am creating an ArcTool and wanting to have my input parameters encompass a drop down menu where files in the TOC can be chosen. The tool is shown below

I currently have a .csv and .shp sitting in my TOC and would like to have dropdown menus present on the "Stats File" and "Spatial Data" parameters so that I can choose them easily, instead of browsing to the file through the folder connection.

I am receiving an "Invalid Drop Item" message.


texts:
I have implemented DBScan in vb.net (forced to...).  I am using a progress bar to track other parts of this tool.  These all are implemented in For loops.  Therefore, incrementing the tracker is not problem.  Set the step to 1 and check for cancel at the conclusion of each loop.

One key item in DBScan is that it pops n elements off the stack as it generates a cluster.  This means that it is not a one iteration to one element mapping.  I tried reassigning the step value with each iteration, but this results in the progress bar recomputing and redrawing occasionally.

What techniques have others used to implement an progress bar outside of a for loop?  I am not worried about total accuracy, but would like a relatively close approximation that can scale across a range of input data sizes and total cluster sizes.


texts: What is the difference between geopackage and shapefiles?

texts:
I'm trying to figure out what Geopackages are and often I read "it's like a shapefile but not quite." Sometimes I also read that it's the equivalent to the ESRI Geodatabase.

What is the actual difference between Geopackages and shapefiles?


texts:
How can I have a dynamic title of the map, which changes depending on the layer chosen? Is this doable? It would save a lot of time and iterations since I should be producing several maps.

To illustrate, I want the title to get the name of the first layer ticked if several layers are ticked. Any hint will be useful.


texts: Calculating floodplain percentage in parcel dataset

texts: I have parcel data for an entire county and a floodplain that lies on top of it. I would like to isolate the 100 year portion of the floodplain and calculate how much of it is within each parcel.


texts: KML Time Stamp With UTC Conversion Not Working

texts:
So I have a KML Placemark with a <when> tag in it that has the dateTime that I want, however when I load it into google earth, the <when> tag disappears.

Here is my code:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.1">
  <Document>
    <Folder targetId="fld1">
      <Placemark>
        <TimeStamp>
          <when>2014-04-15T01:12:06-5:00</when>
        </TimeStamp>
        <Point>
          <coordinates>-82.345144,27.648110,164.787275</coordinates>
        </Point>
      </Placemark>
     </Folder>
  </Document>
</kml>

As soon as I load it into Google Earth, the placemark appears but the When tag does not. Therefore I dont get the timeSlider.

If I simply use the Z notation at the end of the time, then it shows 8:00 rather than 1:00 something...

Any help would be great. I am relatively new to KML.


texts: Remove white gaps which appear when the local tile map is zoom in-out in QGIS?

texts:
I am using QGIS with local XYZ layer. When I zoom in, zoom out and move the map the layers load slowly. First there is white gaps int the edges, then the map loads and the gaps close. On the other hand when I use other platforms like GoogleMaps, when you zoom in or zoom out they blur the map first the they load the map. This way there is a smoother look and you don't see any white gaps while loading the map.

Is there any way that I can implement this to my QGIS application?


texts: ArcMap, Data driven pages, raster layers

texts: In ArcMap I'm trying to display specific raster layers using data driven pages and I got stuck. I have around fifty raster layers and corresponding pages in data driven pages. I only want specific layers on specific pages (one layer per one page is also acceptable). Normaly when I want to achive this with shp's or FC's I use page dafinition but here it's not possible. Any ideas how to do this?


texts: Heads up digitizing - Polygons overlay differently with two different Satellite images

texts:
I have been digitizing building footprints from satellite image using ArcGIS 10.1 . Field data collections is done with the use of four corner point GPS way point. Due to multipath errors, we are using satellite image to complement it for better accuracy.

I have two different satellite images which one of them is a recent one but lower resolution, while the other one is an old one with a better resolution.

I have tried to use the old one with better resolution in my heads up digitizing. It went well .

My problem is once overlayed with recent image, polygons bounce off few meters to the road.

Polygon Coordinates and Projection:
Projected Coordinate System:    WGS_84_UTM_zone_38N
Projection: Transverse_Mercator
false_easting:  500000.00000000
false_northing: 0.00000000
central_meridian:   45.00000000
scale_factor:   0.99960000
latitude_of_origin: 0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree

Old Satelite image Coordinates
Spatial Reference: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Datum: D_WGS_1984

New Satellite image Coordinates
Spatial Reference: WGS_84_UTM_zone_38N
Angular Unit:Degree (0.0174532925199433)
Datum: D_WGS_1984

What could be the problem and how would you go around it?

Thanks