texts:
Anyone ever have to make a fake hill slope? or a fake digital elevation model (DEM)?  What are the best methods for doing so?

For more specifics, suppose I wanted to make a fake hill slope, that was 200m wide by 300m long.  If I pick a stating height of 150m and an end height of 100m.  How should I go about generating this hill slope?  What tools (R, ArcMap, etc...) would be the best method for generating such a hill slope?  Do any tools allow for random variation of the elevation values across the DEM (height and width) while still maintaining the slope value (or range of slopes) I am shooting for?

Asked another way, how do I create a 3D surface that is reminiscent of a real hill slope?

WhiteBoxDev's answer below is pretty close to being the answer I am shooting for.  My only additional  questions relate to randomly altering the elevation of the hill slope (perhaps with a random error?) to make it more realistic - and not perfectly smooth.


texts: What is the "Long-Wavelength Error" reported for SRTM data?

texts:
Assessments of SRTM topographic products (like Rodríguez et al 2005 and 2006) report that height errors showed coherent long-wavelength signatures and describe long-Wavelength error characteristics.

But although I have exhaustively sought what it exactly means I can't figure out to what "long-wavelength" refers to. What is the long and short-wave components of SRTM errors?


texts: Why are OGR (MySQL / PostGIS) connections extremely slow in QGIS 2.0 - 2.8 (when compared to QGIS 1.8)?

texts:
In our project we have been using QGIS 1.8 Lisboa with a MySQL database placed on a server accessed via local network since two years. Everything is working fine and quick. The vector layers can be edited very swiftly, even when the MySQL server is on internet.

Since the change to QGIS 2.0 this is no more the case (the same applies to 2.2, 2.4, 2.6 and 2.8). Even a simple vector layer of 10 features takes a long time to open. It is not a question of rendering. The response time is extremely long even when opening the attribute table values). Naturally, the more features there are, the longer it takes. We were actually forced to continue using 1.8 for editing and we use 2.6 only for printing. It should be noted that we are using 1.8 and 2.6 as parallel installations on the same computers now, and the difference in speed is enormous (c. 1-3 seconds response time for a simple table on QGIS 1.8 and c. 20 seconds on QGIS 2.6). We have computers running Windows 7, Windows 8, Linux Ubuntu 14.04. The operating system has no effect on the speed.

Any ideas?

Switching between two versions is very cumbersome for our workflow (but still much quicker than working with 2.6 only).

Actualization (30th April 2015):
The problem is not related with network. I have made additional tests with localhost, and the results show the same pattern. The old 1.8 version is immediate in its response time, while the 2+ versions take a very long response time. I have however observed an additional, highly remarkable fact - the speed of the QGIS 2.0-2.8 versions is related to the number of tables in the MySQL database!!! If there is only a single table (with 34 simple polygon records), opening its attribute table in QGIS 2.6.1 takes 2-3 seconds. It is not immediate like QGIS 1.8, but it is pretty quick. Now, if I leave this simple polygon layer opened in QGIS 2.6.1 and I add another tables with spatial information into the MySQL database (yet I do not load them in QGIS, they remain only present in the database), the time necessary to open the attribute table of that simple polygon layer will increase (even if the polygon layer is still the only one layer loaded in QGIS)!!! With 4 tables it takes 5 seconds to open the attribute table, with 9 tables it takes 10 seconds, with 81 tables it takes more than 70 seconds to open an attribute table of a single loaded layer. I just do not understand why this happens only by QGIS 2 and higher??? QGIS 1.8 is in no way influenced by the number of MySQL tables. Can be related with GDAL-tools??? I cannot believe we are the only project that has this problem. Can someone replicate this QGIS/MySQL behaviour? Our configuration of MySQL: version 5.6.21, the tables use MyISAM database engine.


texts:
What are the best algorithms to match segments?

I'm trying to match corresponding segments from two map sources, one less accurate but with segment names, and one more accurate without segment names. I want to semi-automatically apply the segment names to the more accurate map.

The requested algorithm has quite a vague description because a "match" is not well defined, and many factors (orientation, relative length, distance) might have different weight in different scenarios; However, I'm looking for a basic knowledge about the general approaches for handling this problem.

Working implementations for open-source environment (PostGIS, shapely, ...) are warmly welcome.

Sample segments : See description below images.


texts: Arc Hydro Batch Subwatershed Delineation does not contain catchment boundaries/falls short?

texts:
When running the Batch Subwatershed Delineation in Arc Hydro for 10.2/10.2.1, many of the output subwatershed areas do not completely enclose the catchments/adjointcatchments they are are representing.  Instead the subwatershed polygon tends to begin about 1/2 raster cell size above the drainagepoint (which followed all necessary batch procedures before application) and in many cases it appears to run out of area from the nearby catchment boundaries and do not cross over the str grid or drainage point.

Why might this be happening and how do I resolve it to encompass the entire catchment/adjoint catchment areas?


texts: Identify duplicate polygons to enable removal after reviewing attributes

texts:
I wish to label overlapping polygons so that I can manually locate each overlap then review the attributes of each overlapping polygon and decide which suite of attributes to maintain before dissolving the polygons to ensure a dataset free of overlap. I know the manual review is time consuming but I cannot see another way of doing this if you don't know which attributes should be maintained before looking at the overlapping polygons.

I wish to be able to calculate the relational area of specific attributes within the dataset without biasing the analysis due to overlap. Ideally I would like a way to add a field that added a code (e.g. number) that would link each polygon with other polygons that it either entirely or partially overlaps with. Then I can determine which of these polygon's attributes to retain.

I have very large datasets with numerous different attributes. I use Arc10 with win7 64 bit processor and ArcInfo licence.  Any ideas and help is greatly appreciated. Thanks


texts: Deleted rasters still appear in geodatabase in ArcGIS Pro?

texts:
I am trying to remove rasters from a geodatabase in ArcGIS Pro 2.2, but the rasters still appear in the catalog and when I run arcpy.ListRasters('*').

However, when I run the delete_management tool, it says that the raster is not present.

Has anyone run into this before, know know how to solve it?


texts:
I read the GDAL's reference for driver implementation and write new driver, but only work in Release mode. So i need debug some parts of program but new driver don't add to GDAL driver's and when to open the file the returned dataset is NULL.

So, question is: How to implement new driver for GDAL in Debug mode?


texts: Add featureset in graphics layer without iterating

texts:
I have a geoprocessing service that returns a large amount of features (featureSet). I want to attach them in a graphics layer, but I don't want to iterate over all them. I tried using the mygraphicLayer.graphics option so as to set the features (array of graphics) to the graphic layer but it does not work.

Have you got any idea how I could do this?

 function addCityFeatureSetToMap(featureSet) {
    console.debug("featureSet    ", featureSet)
    var symbol = new esri.symbol.SimpleMarkerSymbol();
    symbol.setColor(new dojo.Color([0,0,255]));

    //Create graphics layer for cities
    var cityLayer = new esri.layers.GraphicsLayer({
        graphics: featureSet.features
    });
    map.addLayer(cityLayer);


texts: Expand QListWidget in QDockWidget in a QGIS plugin

texts:
I am working on a QGIS plugin. The plugin has a panel with two QListWidget and buttons, labels, etc. I would like to expand horizontally and vertically the lists if panel resized. I set the hsizetype and vsizetype to Expanding for both lists. Unfortunately lists are not expanding vertically after a given size. What should I change in the Qt ui file to expand lists?

You can find my ui file here: https://github.com/zsiki/thematics/blob/master/thematics_dockwidget_base.ui

Additional information

Following Kadir's answer I had to change verstretch to 1 for the main widget  (QDockWidget) to get the desired behavior of my panel.


texts: Displaying MGRS grid when area crosses 2 UTM Zones in ArcGIS Desktop?

texts:
I need to map an area that crosses 2 UTM zones with ArcGIS for Desktop. After having defined the Data Frame projection according one of the two UTM zones involved, I have noticed that the MGRS labels are wrong for the part of the map covered by the other UTM zone.
For example, in this map of central Haiti crossing UTM Zones 18 and 19 I have applied the projection UTM Zone 18N. As you can see, MGRS label are correct for all the 18Q grids, and wrong for the remaining 19Q ones.

Correct MGRS coordinates can be easily checked here: https://mappingsupport.com/p/gmap4.php?tilt=off&mgrs=14SPG34308382&z=5&t=t1

Is it possible to fix this issue, or at least to find a good workaround?


texts: How to change the OpenStreetMap language (or at least the alphabet)

texts:
I just arrived in Belgrade. I am browsing the map in OpenStreeMap.org and in the OSMAnd programme. Everything is in Serbian, street names, points of interest, etc. Check it by yourself.

The difficulty for me with Serbian is that it is coded in a different alphabet, the Cyrillic. This means that searching for a keyword like "museum" returns no results. Even if I knew how to write "boulevard" or "restaurant" in Cyrillic, I would need a different keyboard in my laptop.

Is there any possibiliy of changing the labelling language in OpenStreetMap? Or at least of changing the alphabet, so that the label Београд is displayed as Beograd?


texts:
I'm hoping this is an easy question because I'm a pretty novice web dude:

I've made a visualization in CartoDB, but rather than simply embedding as an iframe, I'd prefer to get the javascript code of the visualization, so I can tweak it once I post it on my blog (to add buttons, etc.).

Where do I get the code to post into my blog page?


texts: Ogr2ogr projection problem- features are shifted 12 miles to the north?

texts:
I have discovered a problem related to the output of an ogr2ogr command line call. My resulting polygon features are 12 miles (.18 degrees) due north of where they should be, but the longitude is correct. The polygon dimensions are also preserved.

The desired output is a shapefile that has been projected using EPSG 4326.  The input comes from a postgres database with postgis where the polygon were digitized in EPSG 900913 (the Google Web Mercator prj).

Below is the command line call that I am using:

ogr2ogr -where " gid = 3 " -f "ESRI Shapefile" E:\Web_Data_Storage\dfr_map_shapes\DFRmap_Tract_3_4326.shp -t_srs EPSG:4326 PG:"user=dfr_webmap dbname=denr host=ipaddress password=xxx port=5433" Tracts2

If I convert the coordinates of a corner of a polygon using gdaltransform I get the correct coordinates.

gdaltransform -s_srs EPSG:900913 -t_srs EPSG:4326 -8853154 4197094.5

Why is this happening and how do I fix it?

Thanks.


texts:
I have web client on Apache and WFS service on Geoserver, both on same computer!
Web client is (Apache) running on port 89, and Geoserver is running on 8181 port!
I do following settings in Apache config file:

<Proxy *>
    Order Deny,Allow
    Deny from all
    Allow from `localhost
</Proxy>

 ProxyPass /geoserver http://localhost:8181/geoserver/web
    ProxyPassReverse /geoserver http://localhost:8181/geoserver/web`
<Location /geoserver>
    Order allow,deny
    Allow from all
</Location>

When I try to run client page I am getting error:

XMLHttpRequest cannot load http://localhost:8181/geoserver/wfs. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://localhost:89' is therefore not allowed access.

Are there any settings that I need to do in Geoserver to enable this?


texts: Reclassifying only part of raster keeping rest unchanged?

texts: I'm trying to reclassify only a certain part of a raster...not the entire raster.  Here's the problem.  I have a LULC raster of the entire US.  There are some issues with the raster that I want to fix, namely some cells that should be water are classified as wetlands.  These are mainly off the coast, not inland.  However, some cells that are classified as wetlands should really be wetlands (mainly inland) and I don't want these to change.  I plan on using a ZIP Code boundary file to determine which cells are considered "off the coast" and which ones are not.  I have access to ArcGIS 10.1 and R with the appropriate raster packages I need so a solution in either software would be fine.  Can anyone help me figure this out?


texts: I use pixel inspector in Arcgis  to find cell values. this is very helpful tool in arcgis. I wanted to know if there is such a tool in Qgis which can help me to find the cell values(Except the identify tool).


texts: Sum line segments in overlapping polygons with Spatial Join

texts: I am using the spatial join tool with counties (polygon) as the target features and gas pipelines (lines) as the join features. I want to see how much pipe is in each county. The problem arises when one pipeline segment overlaps into two polygons. It is then adding the length of the pipeline segment to both counties. I have tried other match options even have their center in but that means the pipeline would have to cross the center of the county.  Without splitting the pipes at county lines is there a more reliable way to do this where it would only put the length into one county, possibly the one it overlaps most?


texts: Leaflet clean checkbox in control layers

texts:
I want to remove a layer from the map and uncheck its control layer checkbox. According to the code bellow the layer is removed from the map, however its checkbox remain checked in control layers. Any hint on how I can clean or uncheck the layer checkbox after removing it from the map?

map.on('layeradd', function (event) {
    if (event.layer == layer1) {
        $('#layer1')[0].textContent = 'MyLayer'
        $('#legend-layer1')[0].setAttribute('style', 'visibility: show')

        if (map.hasLayer(layer2)) {
            map.removeLayer(layer2)
        }


texts: Oracle Enterprise GeoDatabase force logging vs nologging

texts:
Our production Enterprise GeoDatabase running on Oracle has been configured at the database level to force logging.  ArcSDE DBTUNE defines many tables and indexes by default to be NOLOGGING for performance reasons.  Is force logging not a good approach for geodatabases?

Some context : Our production database runs in archivelog mode.  When we do a DB restore of our non-production database from RMAN backups of our production databases, these tables/indexes defined as NOLOGGING in DBTUNE run into block level corruption while the logs are applied after restore. The DBA recommends having force logging turned ON so we don't run into issues while applying the archive logs.   Ref : https://chenguangblog.wordpress.com/tag/nologging/

Ref:
https://resources.arcgis.com/en/help/main/10.1/index.html#/What_is_the_DBTUNE_table/002n0000001n000000/
           https://www.orafaq.com/wiki/Nologging_and_force_logging


texts: Is Coarse / Acquisition PRN pseudo random code value hardcoded into GPS receiver?

texts:
From what I read there are 37 ( 1-32 are for traditional GPS code 33-37 are reserved for other use and ground transmitting ) carefully chosen pseudo random sequences 1023 bits long called Coarse/Acquisition and each GPS satellite has assigned one sequence which is broadcast over and over again.

Is this pseudo-random sequence hard-coded into GPS receivers? Or the receiver can recognise such sequences on fly and store for future timing synchronisation - so then even if new sequences appears they are still recognisable.


texts: This might be a trivial question, but please bear with me. So I am using python and Sentinel images and I want to calculate NDVI, which is  (NIR-RED)/(NIR+RED). So in the case where the final result is NAN, this means that this is a zero division right? i.e NIR=-RED? and if NDVI=0 what does this means? because NDVI ranges from (-1) to (1).


texts: Qgis delimited text layer from an excel file converted to .csv

texts: I'm trying to add a delimited text layer using coordinates. The file has been created in excel and then converted to a csv file in the saving. When I put it into Qgis it doesn't recognise the headings- it puts all the cells as headings and doesn't recognise all the data. none of the headings are just numbers. Attached is a screenshot of the data I'm trying to import. Any ideas as to where I could be going wrong?


texts:
Is it possible to use setView when using an alternative CRS?

var map = L.map('map',{crs: L.CRS.EPSG26915}).setView([4970081,477018], 13);
// Returns: Uncaught TypeError: Cannot call method 'latLngToPoint' of undefined

I assume it's related to using a CRS that isn't lat/long.

Thanks!

To clarify in regards to the first comment...

var map = L.map('map'.setView([44.8846225,-93.2203424], 13); // Works fine
var map = L.map('map',{crs: L.CRS.EPSG26915}).setView([44.8846225,-93.2203424], 13); // Same error as above.


texts:
i have the following code to display dynamically generated addresses on google maps.

var geocoder;
var map;
var bounds = new google.maps.LatLngBounds();

function initialize() {
    map = new google.maps.Map(
    document.getElementById("map"), {
        center: new google.maps.LatLng(10.9027636, 49.8988135),
        maxZoom: 18,
        mapTypeId: google.maps.MapTypeId.ROADMAP
    });
    geocoder = new google.maps.Geocoder();

    for (i = 0; i < locations.length; i++) {

        geocodeAddress(locations, i);
    }
}
google.maps.event.addDomListener(window, "load", initialize);

function geocodeAddress(locations, i) {
    var title = locations[i][0];
    var address = locations[i][1];
    var url = locations[i][2];
    var markerImage = locations[i][3];
    geocoder.geocode({
        'address': locations[i][1]
    },

    function (results, status) {
        if (status == google.maps.GeocoderStatus.OK) {
            var marker = new google.maps.Marker({
                //icon: 'https://maps.google.com/mapfiles/ms/icons/blue.png',
                icon: markerImage,
                map: map,
                position: results[0].geometry.location,
                title: title,
                animation: google.maps.Animation.DROP,
                address: address,
                url: url
            })
            infoWindow(marker, map, title, address, url);
            bounds.extend(marker.getPosition());
            map.fitBounds(bounds);
        } else {
            alert("geocode of " + address + " failed:" + status);
        }
    });
}

function infoWindow(marker, map, title, address, url) {
    google.maps.event.addListener(marker, 'click', function () {
        var html = "

" + title + "

" + address + "

Eintrag anzeigen

";
        iw = new google.maps.InfoWindow({
            content: html,
            maxWidth: 350
        });
        iw.open(map, marker);
    });
}

function createMarker(results) {
    var marker = new google.maps.Marker({
        icon: 'https://maps.google.com/mapfiles/ms/icons/blue.png',
        map: map,
        position: results[0].geometry.location,
        title: title,
        animation: google.maps.Animation.DROP,
        address: address,
        url: url
    })
    bounds.extend(marker.getPosition());
    map.fitBounds(bounds);
    infoWindow(marker, map, title, address, url);
    return marker;
}

The Addresses came over this:

 var locations = <?php echo json_encode($addresses, JSON_PRETTY_PRINT);

Now I have 12 addresses, which are generated dynamically and I get an OVER_QUERY_LIMIT error.
Under 12, there is no error.

Do any of you have a tip how to prevent this?

Thanks a lot!


texts:
Hey I'm building a web application and want to use some GIS services for other planets/moons.

In general, I want to find GIS web resources which apply to bodies apart from Earth.

Anybody have url endpoints to a:

Mars tile service
Titan WMS service
Jupiter's Red Spot place finder
Milky Way routing service
Pluto geocoder

I'm more interested in services on the web, rather than data I would need to download and tile/host myself.

Thanks in advance for all answers!


texts: R Spatial: Trouble with disjointed polygon that has self-intersections

texts:
I have a shapefile which represents the areas of a city which have higher risks of flooding than others. This is a spatial polygon object however the polygon is disjoint and not connected. I also have a database of house locations in the city, and I want to develop a dummy variable to determine whether a house is located in any one of the disjointed parts of the polygon.

I've used the gContains(polygon, points) however I get the error

Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) :
  TopologyException: side location conflict at -9.2199117665562582 38.693948116552484

Which I imagine has to do with the fact that it is not one single polygon but multiple fractions of the same shapefile polygon spread out over the city. Here is a picture of the problem I'm having. The blue areas representing the disjoint polygon and the points being those which I would like to determine if they are located in any of the shaded blue region. For completeness, I've also included the underlying city map.

Maybe the gContains is not the best function to use, but I haven't been able to come across anything to help me out. Happy to share the shapefiles and selection of the points database if that helps.

After running the "Check Geometries" in ArcGIS I see there are several self-intersections of the polygon. I ran the "Repair Geometries" and exported the resulting repaired shapefile to upload into R, however still get the same error message as before. Since these shapefiles were obtained from an external source, it may be that the files are not built in such a way to facilitate the type of analysis I am working on.

I was able to solve the problem with a better understanding of what was occuring with the help of Self intersections for polygons, how to solve?

Exporting the shapefile after running the Union code in ArcGIS removed the errors I was having in R.


texts:
I have a leftleat map and i would like to filter the featuregroup i am displaying.

exp_colJSON = new L.geoJson(exp_col,{
onEachFeature: pop_col,
pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, {
        radius: feature.properties.radius_qgis2leaf,
        fillColor: feature.properties.color_qgis2leaf,
        color: feature.properties.borderColor_qgis2leaf,
        weight: 1,
        opacity: feature.properties.transp_qgis2leaf,
        fillOpacity: feature.properties.transp_qgis2leaf    }).addTo(map).bindLabel(Autolinker.link(String(feature.properties['name'])), { noHide: true });
    }
});

feature_group.addLayer(exp_colJSON);
exp_colJSON.addTo(map);

I also added a combobox and i would like to filter the layer's features by the value selected in the combobox. I want to campare the selected value in the combobox with the "level" property of every feature in the layer.

I have coded so far:

$(document).ready(function(){ // ran when the document is fully loaded
  // retrieve the jQuery wrapped dom object identified by the selector
  var exp_getxocolcolegiosJSON = {};
  var sel = $('#niveles');
  // assign a change listener to it
  sel.change(function(){ //inside the listener
    // retrieve the value of the object firing the event (referenced by this)
    var value = $(this).val();

    //do whatever needed here to filter the layer

    L.geoJson(exp_col, {
    filter: function(feature, layer) {
         return(feature.properties.nivelmodelo.match(/.*value.*\/))
     }
    }).addTo(map);

    // I know this has no sense but i don't know how to code it

  }); // close the change listener
}); // close the ready listener

Regards,


texts:
When using a SQL database as the store, Geoserver requires that the spatial data are stored in a geometry, not geography column. Naturally, as we are dealing with geospatial features and not just Euclidean planar geometry, one would expect the opposite to be true.

Is there any reason for this? Am I missing something important about these two datatypes?


texts: How to update a buffer view with values from another table?

texts:
So basically there's a buffer view defined by:

CREATE OR REPLACE VIEW niceview2 AS
select 'test'::text AS route_type,
st_setsrid(st_union(st_buffer(a.geom, 25::double precision)), 27700) AS geom
from testdata a

I need to have several values counted in it from a few tables. I've had success using:

(select sum(count) from testpoints b where st_contains(st_setsrid(st_union(st_buffer(a.geom, 25::double precision)), 27700), b.geom) ) as units

problem is that this will apparently recalculate the buffers and then do the st_contains. And I need to do it 5 times for 5 different columns so this layer would be really slow

Is there a more efficient way of achieving this outcome?


texts: How to get turn angle, line orientation and line length in ArcMap 9.3.1?

texts:
I am mapping and analysing animal movements using arcgis 9.3.1 along with Hawthorne Beyer's GME, and am trying to utilise the correlated random walk commands, but in order to do this I need to find out the turn angles, step length and line orientation.
Is there a way to easily do this in arcmap for an entire vector layer? Such as is done with calculating total or portion line length by adding a new field and calculating geometry?
I am not adept at using scripts, so if the solution requires this, can you also include instruction on how to run it?
Very many thanks in advance,


texts:
I'm running a .bat script (included) in the cmd prompt. I get the error "ImportError: No module named osgeo" when it tries to call the gdal_merge.py command.
I think this has something to do with the way my paths are set up. I've tried to play around with making sure the .py scripts are in the correct folders. I've read through other similar problems/solutions but mine seems to be slightly different.
The gdal_merge.py script exists within the 'path' directory, however I'm not sure what the connection is between the 'osgeo' module it can't import. Where should the module exist so that it can be imported correctly?


texts: osm2po vs osm2pgrouting: workshop pgRouting

texts:
As osm2po was a more reliable way to convert an .osm I went for that solution to load my date in postgres.

In the workshop however the author points out to check the result after osm2pgrouting.   The following tables must be created:

"classes"
"geometry_columns"
"nodes"
"spatial_ref_sys"
"types"
"vertices_tmp"
"vertices_tmp_id_seq"
"ways"

After the sql (generated from osm2po) has run, only the tables
"geometry_columns" and "spatial_ref_sys" are created.
How do I fix this difference?
Is osm2po incomplete?
Did something go wrong?

Can anyone point me in a better direction?

Thanks alot!


texts: SQL Server table not updating in ArcGIS Online map

texts:
I have a SQL Server table in my ArcGIS Online web app.
It is related to 2 other feature classes also in my map.

The source table is a non spatial table being related on a common field.
Not sure what more to add about the connection.
I've connected to the SQL Server table in ArcMap and related it to my feature classes then published it to ArcGIS Online

The SQL Server table is updated frequently, however these updates to the table are not reflected in my web app.
The updates are new values in existing columns.

I can see them in the .mxd locally, but they do not appear in the SQL Server table being used in my web app. I have not registered the connection on the ArcGIS Server host so the data is copied each time I publish the map service.

What can I do so that my SQL Server table updates with any new values entered into it from Microsoft Access?


texts: Explanation of vector field marker and its options in QGIS 3

texts:
When I go to the symbology of a point feature layer I can change the symbol layer type from a simple marker to a vector field marker. What does the vector field marker do and in what situations is it useful?

Let's say my point feature has a field called "id" and I have 5 points with "id"'s from 1 to 5.

When I add the field "id" as X-attribute I get a horizontal line of markers. The higher the "id" value the longer the horizontal line.

But how long is the line? And what's the purpose of the other options like vector field type, the direction of rotation or the angle unit?

I am looking for good explanation websites. I couldn't find any.


texts: How to create WMS service with GeoServer?

texts:
I am trying to build first WMS service with GeoServer. I have made layer from shape and layer preview works fine. How can I open this layer with QGIS? How to do map with OpenLayers?

I have tried to search GeoServer/OpenLayers guide for Dummies, but not find yet.


texts: Creation of a GIS Competency matrix for employers

texts:
I am in the process of developing a technical GIS competency matrix for new employees. The matrix will be used not only to assess new employees but also to monitor employee development.
The matrix should contain everything from general GIS concepts to advanced server-side GIS and Web-GIS development.

Has anyone had any experience with such a matrix or can shed any light on how this would look?

My first attempt would look this this.

General GIS concepts: CRS models and transformations, Formats, spatial analysis, georeferencing, heads-up digitizing, symbology, colour
Spatial databases: ESRI, Postgresql/PostGIS, Mysql
WebGIS: O-O Programming, Openlayers, Serverside GIS, Gdal, Geoserver/Mapserver/Deegree
General programming: Python, JavaScript, HTML/CSS, JQuery, PHP
Linux Administration: Shells, scripting, commands, administration, monitoring

I am also interested in finding out what kind of weighting these various criteria should carry. My personal bias is that I think I place more value upon the contextual nature of the GIS-Work than the Technical side, i.e that I find it more important to first understand the data before one develops the tools, rather than first to achieve the technical prerequisites and then try to understand what the data or tool actually represents. Due to this, I would rather employ a planner who learnt GIS-Skills, than a programmer who has later developed GIS skills. Obviously this only applies to GIS projects of a strong contextual nature. If I wanted a GIS-Administrator to only monitor JAVA projects on Linux, then i´d obviously prefer an IT professional.


texts: ESRI JavaScript API Web Application and SQL Queries in ASP.NET?

texts:
My workflow (and security restrictions) for a new html/javascript based web application is making me think that I could program an ASP.NET application, inside which I load the 'map' and the 'legends' objects etc using the usual Dojo libraries.

But before I go too far with this approach, I need to know whether usual ASP.NET Sql queries would be enough to kind of feed the map all the information which typically the ESRI Javascript API based applications use. Specifically, currently:

http://myserver/myservice/MapServer/0/query?f=json&where=

could become a standard ASP.NET SQL Server query targeting the database's tables; I suspect I will need to somehow make the return as JSON output and then feed to the web application?

Am I on right track? Has anyone done something like this before? What are the pitfalls I need to be aware of?


texts: Making line thickness scale proportional to zoom in ArcGIS Desktop?

texts:
I'd like to ask how I can see my outline at the same size independent of my zoom level. I've used GeoMedia in the past where you can set the final scale of your map first (let's say 1:5000) and then define the line thickness (let's say 3pt). The thing is that there is an option to draw the line thicker if I zoom in the data view. For example, if I am in the Data View and zoomed to 1:2500, my line will be 6pt thick.

This doesn't affect anything in the final map, but the line will be thicker, which is very nice for looking if features intersect each other and it's easier the line thickness directly instead of going back and forth between the final map scale and the scale I am working on (usually much more zoomed in).

TLDR: How can I configure ArcGis 10 so that my line thickness is proportional to the scale in the Data View?


texts: Use parameter in Field Calculator using QGIS Modeler 3.16

texts:
I would like to set a parameter in QGIS Modeler, which I can use as input in de Field Calculator.
For example:
parameter = GVG (is a string parameter)

Field calculator:
CASE
WHEN Grondwaterdiepte > GVG_min_su THEN 0
WHEN Grondwaterdiepte < GVG_max_su THEN 0
WHEN Grondwaterdiepte < GVG_min_su AND Grondwaterdiepte > GVG_min THEN 1
WHEN Grondwaterdiepte < GVG_max AND Grondwaterdiepte > GVG_max_su THEN 1
WHEN Grondwaterdiepte < GVG_min AND Grondwaterdiepte > GVG_max THEN 2
ELSE 999
END

I would like to replace the GVG's in the Field Calculator by a set parameter so that I don't need to change a lot in de calculator when I am going to use different parameters in the same equation (and only need to change the parameter itself). However, there are some additions to the parameter to assess different attributes.
I tried concatenating the parameter and additions, by many variations like: concat(@parameter, _min_su) or parameter(concat(@parameter, _min_su), but none of them work. How can I make sure that I only need to change this parameter and not all the text in the calculator as well?


texts: Selecting geometry from PostGIS compatible with Open3D

texts:
I have a cube stored as geometry in a PostGIS DB.
I get the coordinates using the following SQL:

select ST_AsEWKT(geom) from polyhedral.vol_1 where id = 1

This presents the following result.

POLYHEDRALSURFACE(
  ((0 1 0,0 0 0,1 0 0,1 1 0,0 1 0)),
  ((1 -0 1,0 0 1,0 1 1,1 1 1,1 -0 1)),
  ((1 0 0,1 -0 1,1 1 1,1 1 0,1 0 0)),
  ((0 0 0,0 0 1,1 -0 1,1 0 0,0 0 0)),
  ((0 1 0,0 1 1,0 0 1,0 0 0,0 1 0)),
  ((0 1 1,0 1 0,1 1 0,1 1 1,0 1 1))
)

would like to use this within the Open3D library.
For .ply files it is simple mesh = read_triangle_mesh(filename.ply)

Any ideas about how I would get the PostGIS output to be compatible with Open3D?

Using the PostGIS python library I can get the result above into python like this.

conn = <connection details>
cur = conn.cursor()
print("connected")
cur.execute('SELECT ST_AsEWKT(geom) FROM polyhedral.vol_1 LIMIT 1')
geom = cur.fetchone()[0]
geom

This parses the query result into a pandas dataframe

s = geom.replace('POLYHEDRALSURFACE', '')
s = s.replace("((","|")
s = s.replace("))","|")
s = s.replace(")","")
s = s.replace("|,|",",")
s = s.replace("|","",1)
s = s.replace("(","")
s = s.replace(s[len(s)-1], '')
s = s.replace(",",";")
s = s.replace(" ",",")

test = s.split(";")

df = pd.DataFrame([sub.split(",") for sub in test], columns=['x','y','z'])


texts: How to Construct Spatio Temporal Clusters of Time Series Data in R

texts: I have count data of one disease for various districts of a state of India on monthly basis. Is it possible to draw spatio temporal clusters in R? I have checked some R- packages for example surveillance and got a paper which states a model 'hhh4' in Section 5 which nearly talks about my problem. But I am not good at data objects of various geographical classes. Is there any way from which I can use my csv file (which contains months, case count and lon lat) and analyse significant clusters over time (in R, secondarily ArcGis)?


texts: How do I allow users to download a GeoTiff on my web map

texts:
I am hoping to get help with what I need to do in order to allow a user on my OpenLayers map to download a portion of a GeoTiff they are viewing through a WMS. The data is hosted through GeoServer. I would like the user to be able to simply click a button and download the bounds of the original file

I am not sure what steps I need to take in order to begin this process. I have worked a lot with serving geospatial data. But never in a way which allows the user to download the data. Any guide in the right direction would be great. Thanks!


texts: What shapes are inside another shapefile

texts:
I am working with some shapefiles and ArcMap. I am new to both of them.

I have a shapefile that contains all of the precincts of Arizona and I also have 9 different shapefiles each of them holds a District. In the real world each district holds multiple districts my current goal is to find out to which district each precinct belongs to.
I learned that ArcMap has a lot of tools that I can use to do the processing of shapefiles. However I am not sure what tool to use or how I can accomplish this.

Each precinct has a unique id and each district also has a unique id.


texts: image rectification with Mercator projection

texts:
I have written an image rectification function which allows a user to select 2 points on an image and 2 points on a map. The function will then rectify the image so it is translated, scaled and rotated to fit the two points on the map.

The code I have written works perfectly along the equator. However this function 'breaks' when moving further north or south. The rotation and the scale are both wrong. I believe this is because of the Mercator projection I use which stretches along the latitude axis.

I am unsure how to proceed. Is there a bit of math I am unaware of which will compensate for the stretching along the latitude? Or am I wrong in my assumption this is because of the Mercator projection?

Below you will find the code I use:

// transformFromPoints is an array with the two coordinates on the image
// They are in lat long
var fromPointA = this.transformFromPoints[0];
var fromPointB = this.transformFromPoints[1];

// transformToPoints is an array with the two coordinates on the map
// They are in lat long
var toPointA = this.transformToPoints[0];
var toPointB = this.transformToPoints[1];

// This is the img we want to transform
var img = this.activeBaseMapImg;

// First update the anchor so the anchor is exactly at fromPointA
// The anchor defines the point around which the image will be drawn
// This is also the point around which the image is rotated.
var cornerPoints = img.getCornerPoints();
var xAnchorPos = M.Math.projectPointOnLine(fromPointA, cornerPoints[0], cornerPoints[1]);
var yAnchorPos = M.Math.projectPointOnLine(fromPointA, cornerPoints[0], cornerPoints[3]);

var xAnchorDir = M.Math.subVector2(xAnchorPos, cornerPoints[0]);
var xLineDir = M.Math.subVector2(cornerPoints[1], cornerPoints[0]);
var yAnchorDir = M.Math.subVector2(yAnchorPos, cornerPoints[0]);
var yLineDir = M.Math.subVector2(cornerPoints[3], cornerPoints[0]);

var xAnchor = xAnchorDir.x / xLineDir.x;
var yAnchor = yAnchorDir.y / yLineDir.y;

// Set the anchor at the fromPointA location
img.setAnchor({x: xAnchor, y: yAnchor});

// Move the image to the toPointA location
img.setPosition(toPointA);

var fromDirection = M.Math.subVector2(fromPointB, fromPointA);
var toDirection = M.Math.subVector2(toPointB, toPointA);

// Determine scale
var scale = toDirection.length() / fromDirection.length();

// Determine rotation
var fromAngle = Math.atan2(fromDirection.y, fromDirection.x);
var toAngle = Math.atan2(toDirection.y, toDirection.x);

var angle = fromAngle - toAngle;

// Perform rotation and scale
img.setRotation(img.rotation + angle);
img.setSize({width: img.size.width * scale, height: img.size.height * scale});

If you need any more information to help in answering this question, please let me know.


texts: Given WGS84 bounding box, is there service that returns most accurate projection to use?

texts:
I'd like to find a service to convert lat/lon in WGS84 to an arbitrary coordinate system in metric that will be the most accurate for my bounding box.  Does something like this exist?

For instance,

lat = {62.0, 61.0} , long = {-149.0, -150.0}

would perhaps bring up Alaska State Plane 4.


texts: Greek characters are not displayed correctly on attribute table and db manager

texts:
I have a shapefile with some churches in Greece and the characters are not displayed correctly on QGIS, although they do displayed correctly on ArcGIS. Also the other shapefiles I have with Greek characters are displayed correctly. I tried to add the shapefile multiple times on ISO-8857-9, Greek, UTF-8 and Windows-1253, but with no result.
Any other way to change that?


texts: Can I use Openlayers Plugin with a different Datum/Projection?

texts:
There is probably already a post for my problem, but after long searching, I could not find a solution for my issue.

I am working with QGIS und I am using the following projection/datum:
MGI / Austria GK East (EPSG:31256) (+proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=0 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs)

I just want to add a layer form the openlayers plugin. For example Google satelite.
the Google satelite uses the following projection/Datum:
WGS 84 / Pseudo Mercator (EPSG:3857) (+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs).

The Problem that I am facing with, is that I can not get the Google layer in the same projection with my other layers (EPSG:31256)! They are not overlapping correct!

Is there a way to do this?


texts: Problem with attribute field in rasterizing tool QGIS

texts:
I am trying to rasterize a vector layer using the "rasterize tool" in QGIS...
The problem is that there is no "attribute fields" available to choose from (even though they exist!)...


texts: Applying rotation before FeatureToRaster_conversion?

texts:
I am converting polygon shapefiles to rasters. With arcpy I am using the function FeatureToRaster_conversion.

gp.FeatureToRaster_conversion(InFeatures, InField, OutRaster, InCellSize)

This function gives control of the cell size but nothing concerning the rotation.

I would like the output raster with a given orientation other than North. Say for example having the picture at 45 degrees.

Any ideas how I can achieve this? With or without using FeatureToRaster_conversion.


texts: I am using ArcGIS server and I have an issue with RuntimeManager.BindLicense(ProductCode.Server) on x64 configuration (dll file name is "ESRI.ArcGIS.Version.dll"), the execution process exits without the possibility to catch an exception. The same BindLicense works perfectly OK on x86 configuration, however this configuration is not suitable because x64 dlls cannot be loaded. My colleagues do not have this issue with x64 configuration, and we could not find any differences yet which could cause this issue. ArcGIS server manager is running fine, I can see it in the web. Any ideas what might be wrong and how to solve this issue? If it is of any help here is the short version of the dumpfile_results.


texts: QGIS how to center the map while drawing a line

texts:
I have a line shapefile which represents a delivery route.
I am able to extent this route by drawing, using an OSM layer for snapping. (works great!)
Since I need a certain zoom level to see what I am doing, I run out of map soon.
Then I have to drag the map manually. And after that go back to drawing.
In the mean while I haven't saved my layer.
Is there a way of centering that map on my last vertex?
While drawing that is?
It wil be in memory somewhere.
QGIS is aware of the last vertex, as shown in my  image.


texts: Delete all fields with the string "RAWI" using QGIS Python

texts: I would like to delete all fields with the string "RAWI". For instance the field "RAWI_WA" and the field "RAWI_BO" should be deleted but not the field "NSG_RAW"


texts: How to develop a web mapping app to upload, process, and download results?

texts:
I would like to make a web mapping application where a user can import their own x-y-z data file, then the application applies filters that I write to it, contours the resulting data, then displays the results and exports to common GIS formats.

I have medium-level web programming skills (PHP, MySQL, etc.).

My question is this: what tools would be the best way to get my end product? Would something like MapScript work? Or do any commercial solutions come to mind?

Thanks for any suggestions anyone might have.


texts: Plot (latitude, longitude) data into a picture

texts:
I have a matrix where en each position I have certain latitude, longitude and each position has a certain value, i.e. amount of rain in that location. I want to plot in C++ the values in an image but right now I am plotting the values based on rows and columns of my matrix, which  does not reasamble the result I would get mapping the value based on (longitude, latitude).

How can I transform between (latitude, longitude) to pixel position in my image?

My objective would be to create some sort of scatter plot where for each transformed (lat,long) I print as a colormap my value.

As I am doing some sort of heatmap I am using a C++ package for that and my current implementation is:

void Data::plotImage(int Nx, int Ny, const std::vector<short int>& image, const std::vector<Spherical> &sph, int counter){
    heatmap_t* hm = heatmap_new(Nx, Ny);
    short int value;
    int aprox_row;
    int aprox_column;
    obtainLimitsMap(sph); // gives the maximum/minimum latitude and longitude as max_lat, min_lat, max_long, min_long
    for(int i = 0 ; i < image.size() ; ++i) {
        value = image[i];
        aprox_row = (int)round(Ny*static_cast<double>(sph[i].latitude-min_lat)/(max_lat-min_lat));
        aprox_column = (int)round(Nx*static_cast<double>(sph[i].longitude-min_lat)/(max_lon-min_lon));
        heatmap_add_weighted_point(hm, aprox_column, aprox_row, value); // TODO: Should plot based on (lat,row) not (row,column)
    }
    unsigned char im[Nx*Ny*4];
    heatmap_render_default_to(hm, &im[0]); // Render map into picture
    heatmap_free(hm); // I do not need the map anymore
    writepng("heatmap"+std::to_string(counter)+".png", Nx, Ny, im);
}

but the transformation I am making is too simple, just considering max/min longitude and latitude and find the corresponding relationship with the actual point.

Moreover, I do not know whether it helps or not but the data I get from the problem is a matrix with data and 2 formulas to convert each pos from the matrix to spherical (phi,lambda). Then I multiply for 180/pi to get the (latitude,longitude) in degrees and now would like to plot it in an image.


texts: How to add a new Parcel inside Parcel Fabric from existing parcel features (ArcObjects)

texts:
I want to create a new parcel in a parcel fabric layer. I have the x, y coordinates for the parcels to be created. I searched around a lot and found really nothing that actually "CREATES" a parcel.

Here is what I have tried so far:

IPointCollection pts = feat.ShapeCopy as IPointCollection; //an existig polygon from another featureclass used for copying x,y,z..
ICadastralPoints cadPoints = new CadastralPacketClass();
for (int i = 0; i < pts.PointCount; i++)
{
    IPoint p = pts.get_Point(i);
    IGSPoint gsPt = new GSPointClass();
    gsPt.X = p.X;
    gsPt.Y = p.Y;
    gsPt.Z = p.Z;
    gsPt.Type = -1; //just tried when nothing worked..made no difference..
    cadPoints.AddPoint(gsPt);
}

IEnumCELines ceLines = new EnumCELinesClass();

ICadastralLines cadLines = new CadastralPacketClass();
ISegmentCollection segments = feat.ShapeCopy as ISegmentCollection; //an existig polygon from another featureclass used for getting segments
for (int i = 0; i < segments.SegmentCount; i++)
{
    ISegment seg = segments.get_Segment(i);

    IGSLine gsLine = new GSLineClass();
    double bearingDeg = utils.GetBearing(seg.FromPoint, seg.ToPoint);
    gsLine.Bearing = (Math.PI * bearingDeg) / 180;
    gsLine.Distance = seg.Length;

    ceLines.Add(gsLine);
}

IEnumGSLines gsLines = ceLines as IEnumGSLines;

IParcelLineFunctions parcelLineFunctions = new ParcelFunctionsClass();
IConstructParcelFunctions constrParcelFunctions = new ParcelFunctionsClass();
ICadastralPacket m_EditPacket = cadExtension.GetEditPacket();

IEnumGSParcels outputParcels = null;
int peRrorLinePointfrom = -1;
int peRrorLinePointTo = -1;

IGSPlan gsPlan = gsParcel.Plan; //gsParcel is an existing parcel..

constrParcelFunctions.ConstructParcelsFromLines(gsPlan, cadPoints, gsLines, m_EditPacket, 1, ref outputParcels, ref peRrorLinePointfrom, ref peRrorLinePointTo);

IGSParcel resultParcel = outputParcels.Next(); //resultParcel is null now..was expecting this to have some value..

Hope anyone would guide me into right direction.


texts: Merging georeferenced images that have overlapping information?

texts: Having trouble to mosaic two overlapping .tif files. I need to merge these two images together into one resampled raster image. Unfortunately I cannot do this with my ability and resource available to me. Notice the image is georeferenced correctly, but it has sections of overlapping data (some with no values). I have tried to use ENVI (georeferenced based mosaic) and ArcGIS (mosaic to new raster) to remove the non data values that overlap good data. What would you do?


texts: Writing to a text file from arcpy.ListFields()

texts:
I'm trying to replicate output to a text file from an attribute table using the write() method to transfer the values to the .txt file. The program writes the values but groups each row of values into braces [] and returns a u"value" format for each value. Additionally blanks are represented by the "None". Finally, the order of the field values are wrong even though I tried to build an order using a criteria list to refine the arcpy.ListFields()

This is the output I'm now getting:

[100.0, u'', 304604, 1298977.349872, 295805.476068, u'113984001', u'113984001', None, None, 4, u'D472B444-7F3C-40E9-85B1-A6A700EBFBEE', u'', 2030, u'BT-501'][100.0, u'', 32204, 1343091.764101, 194014.427182, u'134531001', u'134531001', None, None, 4, u'D18279D5-34D7-46A7-9342-A6AB00A55F48', u'', 3660, u'IS-501']

I would like to get this:

D472B444-7F3C-40E9-85B1-A6A700EBFBEE,113984001,113984001,1298977.349872,295805.476068,304604,,100,BT-501,2030,4,,,
D18279D5-34D7-46A7-9342-A6AB00A55F48,134531001,134531001,1343091.764101,194014.427182,32204,,100,IS-501,3660,4,,,
1ACC967C-241F-4C07-BB7A-A6AC00B632E3,9666001,9666001,1275486.924555,280097.24814,259088,,100,SH-501B,,2,,,A

Here's my code:

inFields = ["UniqueID","Address_ID","Account_ID","X","Y","Ref_ID","Side","Score","Route","Sequence","Service_Day","Demand_2","Demand_3","Week"]

    outputWS = r"C:\RS_Data\Workspace\BATCH_CO_161_imports\gisdb\layers.gdb"

        finalFC = os.path.basename(layer)
        print finalFC
        day = time.strftime("%m%d%Y")
        arcpy.FeatureClassToFeatureClass_conversion(layer, outputWS, finalFC)
        outPath = 'S:\\Shared\\RouteSmart\\CO_161_Cleanscapes\\To_Tower'
        name = os.path.basename(finalFC)
        dataFields = arcpy.ListFields(finalFC)
        file = os.path.join(outPath, ("{0}{1}{2}".format(day,name,".txt")))
        print file
        f = open(file, "w")
        writeFields = []
        for field in dataFields:
            if field.name in inFields:
                writeFields.append(field)
        for row in arcpy.SearchCursor(finalFC):
            fieldVals = [row.getValue(field.name) for field in writeFields]
            f.writelines(str(fieldVals))
        del row
        f.close()

Is there a way I can specify the order in which the values are written, remove the u"value" formatting as well as the []'s around each row, replace "None" with "" and have it write each row as a new line in the text file?

After taking a mix tips from the answers I've gotten it close. The None's are still there but that's the only thing left to resolve. Here's the modified code:

        arcpy.FeatureClassToFeatureClass_conversion(layer, outputWS, finalFC)[0]
        outPath = 'S:\\Shared\\RouteSmart\\CO_161_Cleanscapes\\To_Tower'
        name = os.path.basename(finalFC)

        dataFields = [i.name for i in arcpy.ListFields(finalFC)]

        file = os.path.join(outPath, ("{0}{1}{2}".format(day,name,".txt")))
        print file
        f = open(file, "w")

        writeFields = []

        for field in inFields:
            if field in dataFields:
                writeFields.append(field)

        for row in arcpy.SearchCursor(finalFC):
            # No longer need to reference name property (field.name)
            fieldVals = [row.getValue(field) for field in writeFields]
            # Replace nulls with empty strings
            fieldsVals = ['' if i is None else i for i in fieldVals]
            # Create a comma-separated string of values
            fieldVals = [str(field) for field in fieldVals]
            out_string = ','.join(fieldVals)
            # Write the string--not the list--to the table
            f.writelines(out_string)
        del row
        f.close()

    arcpy.Delete_management(layer)
    arcpy.Delete_management("OldSeq")


texts: Add attribute from underlying polygon based on maximum shared area

texts:
The image below is of vineyard blocks on different soil types (four larger polygons - all different layers).

I would like to create a new attribute in the vineyard layer called soil type and fill it with the value from the four different soil type layer's attribute called description.

In the case where the soil types overlap the vineyard blocks, the attribute soil type should be where more than 50% of the block lies within the soil type. I have looked, but cannot find a answer.

I have tried geomwithin('Alluvial', 'Description') - a function I get in the field calculator after I installed the plugin refFunctions. However, this only works, of course, if the vineyard polygon is completely within the soil type. I have played around with Union, but I cannot see how that will help me.

Can someone please help me? I would prefer a expression, as I have no knowledge of scripting.


texts:
I am using QGIS 3.14.
I'm using an *.ui file for a custom form in the layer properties dialog and everything works fine.
Once the form is open, i change the size of the form with the mouse from the original size to a new size and close it.
When i reopen the form after restarting qgis the form doesn't open in the original size (stored in the *.ui) but with the modified size.
Where does QGis store this modified form size?


texts:
At the moment I am searching for the most effective method of searching an ITable with one or more conditions (like where clauses).

I can perform it manually by iterating through each row and proof manually if the field "xy" has the value "ab". However, because the table has tens of thousands of row this will not be really effective. Does anyone has a effective and elegant solution?

Thanks for your help!


texts: Reading zipped ESRI BIL files with Python

texts:
I asked the question below on StackOverflow, but haven't received any feedback there yet:

I have precipitation data from the PRISM Climate Group which are now offered in .bil format (ESRI BIL, I think)  and I'd like to be able to read these datasets with Python.

I've installed the spectral package, but the open_image() method returns an error:

def ReadBilFile(bil):
    import spectral as sp
    b = sp.open_image(bil)
ReadBilFile(r'G:\truncated\ppt\1950\PRISM_ppt_stable_4kmM2_1950_bil.bil')

IOError: Unable to determine file type or type not supported.

The documentation for spectral clearly says that it supports BIL files, can anyone shed any light on what's happening here? I am also open to using GDAL, which supposedly supports the similar/equivalent ESRI EHdr format, but I can't find any good code snipets to get started.


texts: Opening spatial weights matrix file (SWM) from ArcGIS for Desktop in Excel?

texts:
I successfully generate spatial weights matrix file (SWM) to get the nearest neighbors for polygons shapes in ArcMap 10.2.2 but I am confused how to open this file its a compressed file and when I try t open it I get error to open the compressed file.

I searched on net and found that this file can be converted to table and open in excel but I didn't find the tool that can do this in the ArcMap any suggestion?


texts: Rasterstats throwing error: Cannot open as GDAL raster

texts:
I'm using the rasterstats zonal_stats function (pip installed), and running it on a .vrt file.

Everything runs perfectly in my Ipython notebook, but when I try to run it in my regular python environment I get the following error:

RasterStatsError: Cannot open '.../africa_merged.vrt' as GDAL raster

The error points to the following lines of the rasterstats source code:

rds = gdal.Open(raster, GA_ReadOnly)
if not rds:
    raise RasterStatsError("Cannot open %r as GDAL raster" % raster)

However I when I run this outside of the rasterstats zonal_stats function, it can open the .vrt file, and doesn't throw an error.

Has anyone else had this problem?

UPDATE: this is happening with regular GeoTiff files as well. So the problem isn't because of the .vrt file type.

Some version info:

Python: sys.version_info(major=2, minor=7, micro=9, releaselevel='final', serial=0)
Osgeo:  osgeo.gdal.VersionInfo(): 1110000


texts: Find closest point from origin in QGIS and display measure line in composer

texts:
I have a location point. there are several constraints areas surrounding the point. Is there a way to calculate the closest point of this constraint from my point and to display the distances as a measurement in print composer.

Further to the comments below see the above image.

The red point is the subject site.  The hatched areas are natural heritage constraints areas such as Natura 2000 site and local designated sites.  For the purposes of reports it is necessary to state (for example) that site X is 7Km from the nearest point of the SPA, 9Km from the nearest point of an AONB, 15Km from the nearest point....etc.  It would be ideal if this could be plotted rather than using a layer and snapping if it was possible to use an "autoCAD type function" where by measurements are displayed as an attribute to the line rather than a lable plotted in the print composer.


texts: Problem in displaying a GeoTiff image in JMapFrame

texts:
I'm trying to display a GeoTiff image with no projection but where the geographical extent (bounding box) is known:

reader= format.getReader(new File(sImageName));
GridCoverage2D gridCoverage= reader.getGridCoverage2D(null);
ReferencedEnvelope referencedEnvelope= new ReferencedEnvelope(-106,
-105,55,56, CRS.decode('EPSG:4326', true));
MapContext map= new DefaultMapContext();
map.setCoordinateReferenceSystem(CRS.decode('EPSG:4326', true));
map.setAreaOfInterest(referencedEnvelope, CRS.decode('EPSG:4326', true));
map.addLayer(reader, rasterStyle);
JMapPane pane= new JMapPane();
pane.setMapContext(map);

It does not work (with GeoTools 2.7.5), I get

java.lang.NullPointerException
        at org.geotools.swing.JMapPane.setTransforms(JMapPane.java:1114)
        at org.geotools.swing.JMapPane.doSetDisplayArea(JMapPane.java:607)
        at org.geotools.swing.JMapPane.doSetMapContext(JMapPane.java:505)
        at org.geotools.swing.JMapPane.setMapContext(JMapPane.java:472)

It used to work in an earlier version of GeoTools (2.7.2) with a simpler version of the code:

reader= format.getReader(new File(sImageName));
GridCoverage2D gridCoverage= reader.getGridCoverage2D(null);
MapContext map= new DefaultMapContext();
map.addLayer(reader, rasterStyle);
JMapPane pane= new JMapPane();
pane.setMapContext(map);
map.setCoordinateReferenceSystem(CRS.decode('EPSG:4326', true));


texts: Looking for a sample dataset of international addresses

texts:
I need a test dataset of 100 addresses in each country's local langauge. This is approximately 10,000 addresses. The don't need to be geocoded or resolved to address components. For example,

1 microsoft way seattle wash xxxxx, or
5814 pontiac dr san jose, ca 95123
Почтовый адрес: 125009, Москва, Театральная пл., 1.

are comlete addresses with enough data. Similar ones for addresses in each world langauge would be fine.

Can I get such a non geocoded, random data set for low or no cost? This is to test a google geocoded address cache (as per the API recommendations)


texts: Get layers and layer metadata by geoserver Restful interface

texts:
We are going to access geoserver by RESTful interface, now we meet some problems:

1 List all layers by data store

When using Web administration interface , Geoserver will list all layers(with status of published or not) for a certain data store:

Now how can we get the list by RESTful interface ?

2 Load the feature type metadata

Geoserver will list all the fields when create/edit a layer:

How to load the information by RESTful interface?


texts: Automating Extract by Attribute using ArcGIS ModelBuilder?

texts:
I have a shapefile with several forest data. One of the field is called forest types.

I would create a shapefile for each forest type so Select tool with an easy SQL Expression, such as "forest types" = 'larch' and then "forest types" = 'spruce', but I would automate this operation with ModelBuilder and not changing every time the attribute for selection.

How can I do this?


texts: vb.net arcobjects validate topology ITopologyExtension and esriTEEventHint not declared

texts:
I am trying to use arcobjects with vb.net to validate a topology (already defined). I am caught even when I use the sample provided by ESRI.

The code I use for accessing a topology is:

Dim topoUiD As UID = New UIDClass
topoUiD.Value = "esriEditorExt.TopologyExtension"
Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD),      ITopologyExtension)
Dim topology As ITopology
topology = CType(topologyExt.CurrentTopology, ITopology)

The code I use for selecting errors is as follow:

topologyExt.ClearActiveErrors(esriTEEventHint.esriTENone)
topologyExt.DelayEvents(True)
While Not ((topoError = eErrorFeat.Next) Is Nothing)
    topologyExt.AddActiveError(topoError, esriTEEventHint.esriTENone)
End While
topologyExt.DelayEvents(False)

I get the following errors:

(1) ITopologyExtension is not declared

(2) esriTEEventHint is not declared

I want to check what is missing in the code? Is it the arcObjects library, if yes, what library should be used.

I am a beginner in arcObjects. Please help me on it, thanks.


texts: How to migrate an extension from ArcGIS 9.3.1 to 10?

texts:
I was given an ArcGIS 9.3.1 extension built in Visual Studio (VS) that needs to be updated to work in ArcGIS 10.  When I open the VS project on my computer (that has ArcGIS 10 installed on it) I get a few library reference errors.  I attempted to remove the references then add them back in thinking that would resolve the error, however the errors are still present.  Why are these references not being viewed by VS?  Here are the errors:

Update

I changed the Target .NET Framework from 2.0 to 3.5 and added in the missing file and assemblies via the MSSetup>Primary output from HT.MSExtension(Active).  This removed all the errors except error 4.  Do I need to re-source ESRI.ArcGIS.ADF.dll to ESRI.ArcGIS.ADF.Local.dll?

Update2

I updated ADF to ADF.Local, however I'm still getting error 4.  Any ideas?

  Error 2   Unable to find source file
  'C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MSSetup\stdole.dll'
  for assembly 'stdole.dll', located in
  '[TARGETDIR]' C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MSSetup\MSSetup.vdproj  MSSetup

  Error 1   Unable to find source file
  'C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MSSetup\ESRI.ArcGIS.Editor.dll'
  for assembly 'ESRI.ArcGIS.Editor.dll', located in
  '[TARGETDIR]' C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MSSetup\MSSetup.vdproj  MSSetup

  Error 4   Could not load file or assembly 'file:///C:/Program Files
  (x86)/ArcGIS/DeveloperKit10.0/DotNet/ESRI.ArcGIS.ADF.dll' or one of
  its dependencies. An attempt was made to load a program with an
  incorrect format. Line 195, position
  5.    C:\1_Projects\MS\Source\MS_ArcGIS-Extension\HT.MSExtension\Forms\msManager.resx 195 5   HT.MSExtension

  Error 3   Unable to find source file
  'C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MS.mxd' for file
  MS.mxd', located in '[TARGETDIR]', the file may be absent or
  locked.   C:\1_Projects\MS\Source\MS_ArcGIS-Extension\MSSetup\MSSetup.vdproj  MSSetup

Thank you


texts: Spatial selections using Python in ArcGIS

texts:
I am trying to perform a spatial query using Python in ArcGIS, here is my code:

I'm getting the following error:

  Traceback (most recent call last):   File
  "C:\Users\acarr\Documents_PYTHON
  SCRIPTS\2016SCRIM\02_SPATIAL_QUERY.py", line 19, in
      arcpy.SelectLayerByLocation_management(ExpNod + "/MERGED_NODES.shp","WITHIN_A_DISTANCE",SCRIMNet,SelDist + "
  Meters","NEW_SELECTION")   File "c:\program files
  (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6812, in
  SelectLayerByLocation
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000368: Invalid input data. Failed to execute
  (SelectLayerByLocation).

I'm guessing this is due to the layers not actually being in ArcMap itself. How would I resolve this?


texts: Allow access to PostgreSQL database only through QGIS plugin

texts: I created a Python plugin in QGIS that manage the connection to my PostgreSQL database and loading layers in a specific structure, I've created also some personalized widgets to edit attribute tables so that I can control entered data. My problem is that I'm afraid that someone access my database with the username and password provided through other software such as ArcGIS, so how can I make access to my database only allowed through my plugin ?


texts: How to configure the Web App Builder Coordinate widget to show Z as well as X and Y?

texts:
Is there a way to configure / edit the default coordinate widget to show elevation (Z) along the X, Y values just as is possible in a 3D application but in a 2D application? If possible, this would be the ideal solution.

See capture 1 showing default OOTB coordinate widget in a 2D web appp.

See capture 2 for closeup of the widget with a coordinate placed on map.

Alternatively, can I obtain XYZ values from a DEM service? We've tried this and although we can get the elevation value from in a popup we don't seem to be able to get the cell X and Y coordinate as well as the elevation cell value in one popup.

See capture 3 for a popup showing an elevation (Z) value obtained from a DEM exported to server as image service.

Having to obtain XY values using a widget then getting the Z value from a popup is a bit cumbersome so I am basically looking for a solution that will allow me to obtain all 3 values (X,Y,Z) from a single click on a map.

Capture 1:

Capture 2:

Capture 3:


texts: A few PgRouting database design questions

texts:
Does someone know how the cost field in PgRouting is calculated?
What are 'source' and 'target' good for? 'osm_source_id' and 'osm_target_id' contain the same information using different numbers.
Is there a good reason to use 'cost' and 'r_cost'? Every cost in my table is smaller than 1. So, 0.00045... could mean both directions and 1.00045 oneway.
When it comes to routing what are 'clazz', 'flags' and the line string good for?


texts: Getting started with QGIS, but can't find sample data files on my machine

texts:
I'm just getting started with QGIS. After installation, the User Manual directs you to load raster and vector data. It says "browse to folder qgis_sample_data/raster/." Well, that folder does not exist under Program Files/QGIS Brighton. It next directs you to find and load an .img file. There are no .img files in the current installation of QGIS.
So, what's going on here? Is it just a bad user manual?


texts:
The shapefiles I was sent did not contain a .prj file, and instead came with these .xml files that seem to be the equivalent of the projection file. QGis does not seem to 'read' these XML files when I import the associated shapefiles, and results in the imported polygons being totally out of place and out of scale.

<?xml version="1.0" encoding="UTF-8"?>
<data>
    <coordinateSystem>
        <name>Universal Transverse Mercator - Zone 51 (N)</name>
        <datum>World Geodetic 1984 (WGS84)</datum>
        <system>Transverse Mercator</system>
        <unit>Meter</unit>
        <majorAxis>6.378137000000000000000000e+006</majorAxis>
        <eccentricity>8.181919084262148600000000e-002</eccentricity>
        <centerX>0.000000000000000000000000e+000</centerX>
        <centerY>0.000000000000000000000000e+000</centerY>
        <centerZ>0.000000000000000000000000e+000</centerZ>
        <rotationX>0.000000000000000000000000e+000</rotationX>
        <rotationY>0.000000000000000000000000e+000</rotationY>
        <rotationZ>0.000000000000000000000000e+000</rotationZ>
        <scaleAdjustment>0.000000000000000000000000e+000</scaleAdjustment>
        <scaleX>9.996000000000000400000000e-001</scaleX>
        <localScaleX>1.000000000000000000000000e+000</localScaleX>
        <scaleY>9.996000000000000400000000e-001</scaleY>
        <localScaleY>1.000000000000000000000000e+000</localScaleY>
        <falseEasting>5.000000000000000000000000e+005</falseEasting>
        <localOffsetX>0.000000000000000000000000e+000</localOffsetX>
        <falseNorthing>0.000000000000000000000000e+000</falseNorthing>
        <localOffsetY>0.000000000000000000000000e+000</localOffsetY>
        <centerLat>0.000000000000000000000000e+000</centerLat>
        <centerLon>1.230000000000000000000000e+002</centerLon>
    </coordinateSystem>
</data>


texts:
A python script was written about 18 months ago by a person who has now left.  It produced the required outputs then.  I've been asked to run it again but with different (finer resolution) data inputs.  The input dataset has been split into 20 sub-sets of approx 2,700 data points each.  However, the script crashes ("python.exe has stopped working") after approx 300 data points have been processed (range 295 to 306 and does NOT always fail on the same record).

As its old(ish), the script was written using arcgisscripting and not arcpy. Broadly it does the following using cursors:

For a given point, calculate the cost distance (using gp.CostDistance_sa) with a cutoff of 60 minutes travel time.
Calls gp.ExtractValuesToPoints_sa to extract all the individual values at each data point and outputs a feature class to a file geodatabase.
Reads the feature class created in b) above and writes the values to a CSV file (omitting any points with "No Data" (value -9999)).

Repeats 1, 2 and 3 for all remaining data points in the input file.

Processing time is approx. 1 minute per data point on average. Here are some relevant technical specifications:

The PC has a quad core  Intel i7-2720QM CPU running at 2.20GHz with 8GB RAM running Windows 7 (64 bit).
Python version is 2.6.6 (shell also states "[MSC v,1500 32 bit (Intel)] on win32).
ArcMap 10.0 (SP4) is also installed.

I've tried running it on a different PC (so far without crashing).
Currently the job is running successfully (but more slowly) on an older PC and has reached 419 records without crashing.  The relevant specifications for this machine are:

Intel Core 2 DUO E7500 processor running at 2.93GHz with 4 GB RAM and 64bit Windows 7.
Python version 2.5.1 (shell also states "[MSC v,1310 32 bit (Intel)] on win32).
ArcMap 9.3 is installed (no mention of any Service Packs).

Can someone offer some advice about why the script seems to work for a while then crash and how to resolve it?

The fact that a different PC appears (so far) to handle the script suggests something "environmental".

As an update, the PC running ARCGIS 9.3 is still successfully processing the data and has reached 1,300 data points processed (and still counting).  A colleague also ran the data on their PC running ARCGIS 10.1 - it crashed after 267 records on two separate occasions.  Although not conclusive, the common thread seems to be that Arc 9.3 will process the data but Arc 10.x will not.


texts: How do I create/write a Shapefile with UTF-8 encoded attributes with ogr2ogr?

texts:
I'm using GDAL/OGR's ogr2ogr command line tool to export data from a PostGIS-enabled PostgreSQL database to various GIS file formats, Shapefile amongst them. To create a shapefile with the default encoding (ISO8859-1 a.k.a. latin1, see Which character encoding is used by the DBF file in shapefiles?), I'm using a command line like this:

ogr2ogr \
    -f 'Esri Shapefile' \
    $OUTPUT_PATH \
    -t_srs $OUT_SRS \
    "PG:dbname=${DB_NAME} user=${DB_USERNAME} password=${DB_PASSWORD} schemas=${SCHEMA_NAME}"

But in my data, there may be features with arbitrary languages and scripts in the attribute values. I'd like to preserve these values and thus think that the *.dbf file of the exported Shapefile should be UTF-8 encoded. (The attribute names are guaranteed to be within 7-bit ASCII.)

How can I get ogr2ogr to write a UTF-8 encoded *.dbf file when exporting to Shapefile? The documentation of the GDAL/OGR "ESRI Shapefile / DBF" driver is explicitly ambigous (sic!) about the ENCODING option in the "Layer Creation Options":

  The default value is "LDID/87". It is not clear what other values may be appropriate.

And will the used encoding be indicated in the *.dbf file itself or in an accompanying *.cpg file? (I guess the latter, as UTF-8 probably isn't per se a valid (dBASE DBMS) DBF encoding.) If the latter, will ogr2ogr create the *.cpg file or do I have to create it manually?


texts: Add WFS layer (JSON) in Leaflet using Leaflet layer tree control

texts:
I'm using the layer tree control to manage the display of layers. When I add WMS services, it works fine. But when I add WFS, this is not loaded.

I'm using the Leaflet layer tree plugin code and GeoServer.

Below is the part of the code I use to add WFS layer on my main page:

...{
 "code": "space",
 "name": "space 1",
 "active": true,
 "selectedByDefault": false,
 "openByDefault": null,
 "childLayers": [],
 "selectType": null,
 "serviceType":"WFS",
 "params": {
            "request": "GetFeature",
            "service": "WFS",
            "styles":'',
            "outputformat": "application/json",
            "typeName": "rbrava:space_1",
            "url": "http://10.11.15.170:8080/geoserver/rbrava/wfs"
                      }
                    ...

I used the Chrome inspector and when I click on the layer checkbox to load that it gives me the following error:

"XMLHttpRequest cannot load...
No 'Access-Control-Allow-Origin' header is present on the requested resource.
Origin 'null' is therefore not allowed access".


texts:
I am currently working with SCP within QGIS and am trying to classify the contents of a field (different flowers, grasses, bare ground etc.). I have loaded a raster image of my field and am now attempting to create my ROIS. However, when I try and create a new training file, I am asked multiple times to set the CRS for the field (I want to use WGS 84). Additionally, when I attempt to create a new ROI vertex, I am asked multiple times to set the CRS for the layer and each point. Lastly, I am not able to save the ROI and am prompted with a message saying that I do not have a raster layer loaded and that it cannot be projected.

I believe that I have set up my multiband image and RGB list correctly.

Does anybody know what might be wrong?


texts: Geocode points with append data ArcGIS Online

texts:
I am trying to add points with address data in a .csv file to a hosted layer in ArcGIS Online with the hops of ArcGIS geocoding the points for me. When I try to append the data, I am prompted to provide location information (longitude and latitude) I read this post on the ArcGIS site and it suggests it's possible. https://developers.arcgis.com/rest/services-reference/append-feature-service-layer-.htm
I am able to have points geocoded when I add a new layer, but I would prefer not to have to add a new layer each time I have new data.


texts: Retrieving a Geocoding Result from ArcGIS Online World Geocoder

texts: I ran the World geocoder service through AGOL on a lot of points (I forgot how expensive it is - dumb dumb mistake, alas!). Worse still, the ArcMap Desktop mxd where I had the results open as I was going through the Tied records and manually matching them, died in the middle of the process, before I saved the mxd. Ack!!! I am trying to figure out if I can get those geocoder results back into Desktop somehow, but I don’t know where to look for them. They better be retrievable, given how much they cost me! I have a standalone Standard ArcGIS Desktop license, and an Organization ArcGIS Online account.


texts: Is there tutorial available for PostGIS setup in ArcSDE 10.1 geodatabase on Linux?

texts:
Is there any tutorial for Set up PostGIS in an ArcSDE 10.1 geodatabase on Linux?

There is one in tutorial arcgis 9 for ArcGIS version 9, but I couldn't find any for ArcGIS 10.1.


texts: arcpy.Clip_analysis() Automatically Affixing .shp.shp to Target Feature Class Name?

texts:
I created a Python script that was previously working fine, but is now giving me an error when I try to use the arcpy.Clip_analysis() function in arcmaps 10.2.1.  I am not adding a filetype extnsion to any of my feature class names when creating them.  Previously, this worked fine.  Now however i looks like ".shp" is being affixed to my feature class name twice, making the name invalid.  Here is my code:

    in_features = os.path.join(arcMapsAuto.pc_gdb, arcMapsAuto.merged_feature_class )
    clip_features = os.path.join(arcMapsAuto.shapes_gdb, "Terminal_Polygon_Merge")
    out_feature_class = os.path.join(arcMapsAuto.tc_gdb, "All_Terminal_Clip")

    # Check to see if feature class already exists
    if arcpy.Exists(out_feature_class):
        print "All_Terminal_Clip creation: FAIL."
        print " --All_Terminal_Clip already exists."
        #return

    print in_features
    print clip_features
    print out_feature_class
            # printed out:  {geodatabas}\All_Terminal_Clip

    # Execute Clip
    try:
        arcpy.Clip_analysis(in_features, clip_features, out_feature_class)
        print " All_Terminal_Clip creation: SUCCESSS."

    except Exception as err:
        arcpy.AddError(err)
        print "All_Terminal_Clip creation: FAIL."
        print err
        # printed error:  ERROR 000210: Cannot create output {geodatabase}\All_Terminal_Clip.shp.shp ERROR 000354: The name contains invalid characters

        return

I am missing something?  Is arcpy.Clip_analysis() affixing .shp.shp to my target feature class name, and is this a new issue?  I swear I've used the same code before without a problem?  Any ideas on how I can fix this?

  {geodatabase}\All_Terminal_Clip

  All_Terminal_Clip creation: FAIL.

  ERROR 000210: Cannot create output
  {geodatabase}\All_Terminal_Clip.shp.shp ERROR 000354: The name
  contains invalid characters


texts: Identifying polygon/point shapefiles in PyQGIS?

texts:
I have a script which performs various functions on shapefiles, most of which are polygon-based. However, some shapefiles are point-based and I have to specify where these files are by editing the script. Is there a technique in python which allows the script to identify which shapefiles are polygons, which are points, and run the necessary functions. I have provided a snippet (tested) of my script below which focuses primarily on polygon shapefiles and I have commented out the last line of code which includes the function I use for point shapefiles.

I am using QGIS 2.6.0.

##Test=name
import os
import glob
from os.path import expanduser
home = expanduser("~")

#   Folder path of the directory and the "Results" folder
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

def run():
#   Set directory and search for all .shp files
    os.chdir(path_dir + "Shapefiles\\")
    for fname in glob.glob("*.shp"):
#   Clip .shp files with the polygon shapefile "Grid.shp" and save files to Result folder (How to specify this for polygons only?)
        polygon_output=processing.runalg("qgis:clip", path_dir + "Grid.shp", fname, path_res  + "/"+ fname)
#   (How to identify point shapefiles and run the following command?)
#        point_output=processing.runalg("qgis:distancetonearesthub", path_dir + "Grid.shp", fname, 'HubName', 0, 0, path_res  + "/"+ fname)
run()

I managed to get a working script (even though the entire code could be vastly improved!).

##Test=name
import os
import glob
from qgis.core import *
from os.path import expanduser
home = expanduser("~")

#   Folder path of the Results for shapefiles
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

def run():
#   Set directory and search for all polygon .shp files
    os.chdir(path_dir + "Shapefiles\\")
    for fname in glob.glob("*.shp"):
        shapefile = QgsVectorLayer( os.chdir(path_dir + "Shapefiles\\"), fname, "ogr" )
        print shapefile.geometryType()

#   Clip .shp files with the polygon shapefile "Grid.shp" and save files to Result folder
        if shapefile.geometryType() == 3:
            polygon_output=processing.runalg("qgis:clip", path_dir + "Grid.shp", fname, path_res  + "/"+ fname)

#   Use Distance to nearest hub function on point.shp files with the polygon shapefile "Grid.shp" as Hub layer and save files to Result folder
        if shapefile.geometryType() == 3:
            point_output=processing.runalg("qgis:distancetonearesthub", path_dir + "Grid.shp", fname, 'HubName', 0, 0, path_res  + "/"+ fname)
run()


texts: Is it possible to use JTS for calculations on data in geographic coordinates?

texts:
I have some spatial data that has all of its coordinates as lat/lon pairs (with about 10 digits decimal precision), it's stored in a database as WGS84 data. Some of the data is represented as polygons which are the resulting union of some smaller polygons whose boundaries are stored. I also have a number of points from which I build line segments (just 2 points in each segment) which I use later for intersection tests with the polygons.

I'm using a SpatialIndex to improve my queries so I insert the envelopes of all polygons in a tree (tested with both QuadTree and STRtree). Then, I connect two points into a line segment and I'm using its envelope to query the tree for possible intersections. The problem is that I get pretty much all the polygons as a result which is clearly wrong. To give you some idea about the real scale of my data, I have about 100 polygons that cover the whole North America, each line covers a very very small part of a single polygon. Ideally, I would expect no more than 2 polygons as a result.

I'm using JTS to do this calculation and I'm aware that it's not really suited for spherical data so can you suggest me another library/tool to achieve the desired behaviour or possible a workaround (for example, projecting before using JTS)?

EDIT: I'm using postgresql + postgis extension, storing the all spatial data as geometry in 4326.


texts: Can you recommend a free service for validating a global address?

texts:
I need an API or a database that is freely available. It should allow being included in a commercial application.

Something like Google Places API, but with better licensing terms.


texts:
In QGIS (v3.2) is there a spatial sort similar to ArcGIS (http://desktop.arcgis.com/en/arcmap/latest/tools/data-management-toolbox/how-sort-works.htm#ESRI_SECTION1_9C199CD3D68A41029943D4FF79500A73) where you can sort the spatial attribute as LL, UR etc

The reason for this is that I am creating points and start/end points of roads and if there are segments of the road that aren't in the order of the flow of the road sometimes the xat(0) etc get placed in the middle of the road.

I see a solution from last year in Sort line features into route order in the attribute table and rather than coding this I am wondering if there is a tool now that can be used in modeler. See my other question on this Sort a table by multiple columns in QGIS where the last answer introduces coalesce


texts: Using PyShp to create polygon shapefiles?

texts:
I have a script for creating points shapefiles, but I can not fix it for make a new polyline shp file.

I have obtained the data from a excel file and extract data in three differents lists
Lon,lat,label

import shapefile as shp
from openpyxl import load_workbook
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
r = 'myfile.xlsm'

wb = load_workbook(filename= r, data_only=True, use_iterators = True)
ws = wb['Export_Area']
cell_range = ws ['B5:D9']
tuple(ws.iter_rows('B5:D9'))
lista = []
for row in ws.iter_rows('B5:D9'):
    for cell in row:
        lista.append(cell.value)
        Longitud_area = (lista[1:len(lista):3])
        latitud_area = (lista[0:len(lista):3])
        puntos_area = (lista[2:len(lista):3])

creating the shp file

out_file = 'ExportArea.shp'

x,y,id_no=Longitud_area,latitud_area,puntos_area
logic = [True,False,True]
w = shp.Writer(shp.POINT)
w.autoBalance = 1 #ensures gemoetry and attributes match
w.field('X','F',10,5)
w.field('Y','F',10,5) #float - needed for coordinates
w.field('etiqueta','D') #date

for j,k in enumerate(x):
    w.point(k,y[j]) #write the geometry
    w.record(k,y[j],puntos_area[j]) #write the attributes

creating the prj file.

prj = open(out_file +'.prj', 'w')
proyeccion      ="GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
prj.write(proyeccion)

It works fine for a point shp file but I can not see how to do the loop for making a poly shp file


texts: Transforming Esri JSON to Shapefile without using ArcPy?

texts:
I have a large ESRI json (which I got from a request at an ESRI Rest API) that I need to convert to some sort of vector file, preferably a Shapefile.

It seems there are ways of doing this using ArcPy, but I do not own an ArcGIS license. Is there a way to work around this or is ArcGIS the only way to do it?


texts:
Following the instructions, I've updated the iwave.h and iwave.cpp files using the official spectral responses from DigitalGlobe (In the iwave.h file I don't know if it's necessary to use more accurate wavelength ranges). Now I don't know how (and I don't have much time!) to go on the next steps.. (update geomcond files and recompile grass), I have never compiled a software!! :)
Could someone help me? I can send you all the updated files, of course..


texts: Converting lat/long to EPSG:27700 in geometry generator expression

texts:
I have incorporated a geometry generator into a style in QGIS 3.4 which creates a callout line from a feature to its relevant spatial position.

For reasons out of my control our projects use CRS: EPSG:27700 but some of the layers we need to import use CRS: EPSG:4326. QGIS handles this fine and can plot the features, but the geometry generator cannot.

I know I could convert beforehand but this is not practical for all our users so I just want the expression to handle this for them.

I am looking for a method to convert from '4326' to '27700' within the expression but I am not currently aware if this is possible.

My current geometry generator expression

make_line(make_point($x, $y), make_point("Longitude", "Latitude"))

The "Longitude" and "Latitude" fields relate to the fields in the relevant layer which is projected with EPSG:4326.


texts: Convert measurement geometry from ground reference frame to satellite reference frame

texts:
I need to convert solar_zenith_angle, solar_azimuth_angle, platform_zenith_angle, and platform_azimuth_angle from a reference frame at the WGS84 surface to a reference frame at the satellite observer, in 817km altitude.

Using lat/lon of the ground reference frame, and using the platform_azimuth_angle and platform_zenith_angle, the first step would be to calculate lat/lon coordinates of the satellite ground point.

Next, the four angles can be converted from the ground-based reference frame to the satellite-based reference frame.

Question (a): Is there am existing script/library to do these kind of conversions, preferably in Python?

Question (b): If there isn't, how do I consider the WGS84 ellipsoid in my trigonometric calculations (instead of simply assuming a spherical Earth)?

EDIT 1: Why am I doing this? I want to use a radiative transfer model to simulate spectral measurements from the S-5P satellite. The scenarios I want to simulate
are defined by geometry at the surface. However, the radiative
transfer model I'm using expects the geometry at the observation
platform.


texts: Reordering polygon shapefile permanently?

texts:
I've written a script which produces a polygon shapefile of the extents of all the thousands of rasters in our library. I want to organise the FID so that the files with the largest extents are on the bottom so that an area of interest can be easily selected.

By organise FID - i mean that the layers should be drawn (in the shapefile) from the smallest extent to the largest extent as otherwise (like in layer ordering) the record of the dataset that has the highest FID will get selected over anything below it. This is not a simple (temporary within map/layer) sort but it has to be to permanently change the order in which the rows are written into the shapefile.


texts:
Taking into account topology and existing building heights, would it be possible to do a MCDA and examine where densifying through new construction would be most profitable*, assuming that new construction will utilize the maximum permitted building height?

Profitable defines as to where we would get the highest number of new homes, or building volume in cubic meters if you wish. For simplicity sake: lets treat current buildings like overlays on the ground and pretend new construction can take place on top of them as well.

I've been trying for a couple of days now and I'm starting to think that it's close to impossible with Qgis?

1. Is this kind of advanced MCDA realistic to achive with free software? If so, could anyone outline a simple workflow?

2. Would you go for a raster MCDA or try to use the 50m vector points and solve this vector based instead?

What I have at hand:

Elevation_2m - raster
Elevation_50m XYZ - vector points
3D Buildings - MultipolygonZ-type vector with building heights as z-values
Building height restrictions - vector polygons
Qgis 2.18 Las Palmas, PostGIS up and running with PostgreSQL 9.3.1

And when projected in Qgis2threejs it's seems to be working out. The houses stand on the 2m elevation raster with correct building heights. And the height limitation (an airport approch zone) is the purple "sky" you're seeing:

More details:

Elevation_2m - raster

Elevation_50m XYZ - vector points

3D Buildings - MultipolygonZ-type vector

Building height restrictions - Airport approach zone which restricts building heights. I haven't digitized it completely yet, hence I'm attaching the approach zone that I'm currently digitizing.


texts: Converting layer with symbology using multiple attributes to KML file using ArcGIS Desktop?

texts:
I'm trying to convert a point feature layer with a symbology that uses multiple attributes (a color-scheme attributed to the percentage of positive malaria cases and a symbol size attributed to the number of cases evaluated) into a KML file using the coversion toolbox of ArcGIS 10.1.

Unfornately, when attempting this using the Layer to KML tool, the resulting KML file does not retain the full color-scale defined in the original point feature layer; it represents the symbol sizes correctly, but all the points are reprensented in only 2 types colors (dark and light green). However, a full scale of dark green to dark red should be represented in this KML file similarly to the original point feature layer displayed in ArcGIS 10.1.

I've also tried using the Map to KML tool, but I obtained the same results.

Is it possible to convert a point feature layer with a symbology using multiple attributes to a KML file via the Conversion tool box of ArcGIS 10.1, without losing any elements of the symbology?


texts: How to remove shapefile polygon vertex marker points in QGIS 1.8.0

texts:
I have created a shapefile layer and drawn a polygon to show the boundary of a woodland. How can I edit the polygon to remove the crosses that denote the intersection of two points on the boundary?

You can see that I have created a simple single polygon around the woodland area. I wish to get rid of the red crosses shown. How do I do that?


texts: In QGIS, how to get an angle of a polygon / get longest line in polygon and angle of that line?

texts:
I am using QGIS to output data to a game engine, I use polygon centroids to find center points of polygons of houses to output to the game engine for house placement - now I'd like to take it one step further and get a correct(-ish) orientation of the house.

So my basic idea is to either straight up use the polygons to somehow get the angle relative to true north - or just find the longest side of a polygon house shape and then determine the azimuth/angle of that line - and preferably add this data as a new field to my polygon centroids.

I am still quite new to QGIS and GIS data in general but have been able to make a lot of use out of the GIS data I have.


texts: arcpy.env.overwriteOutput = True not working for one dataset

texts:
I am doing a basic copy from sde database to file geodatabase. I have a loop going through each feature dataset in the sde database, this loop is really basic it just uses arcpy.Copy_Management to copy everything over. I have arcpy.env.overwriteOutput = True at the beginning of the script after importing in arcpy. It does not throw me that error for any other feature dataset, even though they are also overwriting a FeatureDataset with the same name.

the error it throws is  :

  ERROR 000605: Output Sidewalk already exists within
  C:\Users\stack\Desktop\Testing\scripts\Dump.gdb\Sidewalk. Failed to
  execute (Copy).

I even added arcpy.env.overwriteOutput = True to the beginning to of the function. I also put a line in the exception the prints out the status of env.overwriteOutput and it is true at that time when it breaks.

This does not make any sense to me. Has anyone else had this problem?


texts: OGR says that my data source is invalid when I use the method QgsVectorLayer in PyQGIS

texts:
I am begining with PyQGIS and QGIS in general. I wanted a way to create a layer via python, so i downloaded "Plugin Builder" and starded coding. Here's the code :

def activatePlugin (self):
    global limitAer,limitSou ,limitOra , limitENE, limitCha, layer, QgsVectorLayer

    layer = QgsVectorLayer('/Users/fse/Desktop/54_028_380_EXE_CCE_SID_CIR_001_B2/SHAPE/couche_test.shp', 'couche_test', 'ogr')

    QgsMapLayerRegistry.instance().addMapLayer(layer)

    if not layer.isValid():
        print "Layer was not loaded!"
    else:
        print "Layer was loaded successfully!"

I get a "Layer was not loaded" response and a print in the OGR console telling me that the source data is invalid.

I have tried all sorts of slashes, backslashes etc... It is basically always the same.


texts: QGIS clipping raster with CUTLINE_ALL_TOUCHED

texts:
How should the CUTLINE_ALL_TOUCHED option in gdalwarp (referenced here) be used to correctly include all of the touching pixels of a vector boundary to a clipped image?

As can be seen below, the output file has nan and 0 values rather than containing a clipped raster as expected:


texts: Geometry from multiple selected features into one geometry

texts:
Lets see if I can explain what I'm trying to do. First of all I have developed an AddIn for ArcMap that can take one selected feature (Polygon) on a layer, grab the geometry of that feature, and save it to a SQL Server database. What I would like to do is to save multiple selected features, not as individual geometries (ie. one row per feature/geometry) but as one "piece of geometry" containing the selected polygons. This is the code I have so far:

IFeatureLayer featureLayer = ArcMap.Document.CurrentContentsView.SelectedItem as IFeatureLayer;
        if (featureLayer != null)
        {
            IFeatureSelection featSel = featureLayer as IFeatureSelection;
            IEnumIDs idList = featSel.SelectionSet.IDs;

            int index = idList.Next();
            List<int> indexes = new List<int>();

            while (index != -1)
            {
                indexes.Add(index);
                index = idList.Next();
            }

            IFeatureClass featureClass = featureLayer.FeatureClass;
            IFeature feature = featureClass.GetFeature(indexes[0]);
            IGeometry geometry = feature.Shape as IGeometry;
            // Save to database
        }

As you can see I can get all the indexes of the selected features of the selected layer and get the geometry, the challenge is to "concatenate" all the geometries into one...hope it makes sense :)

Any suggestions?

P.S ...and if any of you have a better way of looping through the idList...please let me know :)

UPDATE:

A BIG thanks to Petr! Followed your instructions and it worked in the first try!

Here's the code I ended up with:

IFeatureSelection featSel = featureLayer as IFeatureSelection;

                if (featSel.SelectionSet.Count > 0)
                {
                    ITopologicalOperator resultPolygon = new Polygon() as ITopologicalOperator;
                    IGeometryCollection geometriesToUnion = new GeometryBag() as IGeometryCollection;

                    ICursor cursor;
                    featSel.SelectionSet.Search(null, false, out cursor);

                    IFeatureCursor featureCursor = cursor as IFeatureCursor;
                    IFeature feature;
                    while ((feature = featureCursor.NextFeature()) != null)
                    {
                        geometriesToUnion.AddGeometry(feature.Shape as IGeometry);
                    }

                    resultPolygon.ConstructUnion(geometriesToUnion as IEnumGeometry);

                    // Save resultPolygon to a database
                }