texts: Creating a version to host only proposed infrastructure, that is reconciled and posted when the work is completed?

texts:
I have come across a scenario where I have been asked to find a way to create a version that stems off of a DEFAULT version of ArcGIS Server and hosts only proposed infrastructure data.  This is being done so that there is a way to visualize the proposed infrastructure with the current infrastructure and so we can be ahead of the game when it comes to updating the actual constructed layers.  I am trying to determine whether or not this is possible and if it is, how can this be done without impacting any other data hosted on the server?

I have come Phased Project versions, but I am unsure as to how the work flow would be implemented and whether or not it would work for what I am trying to do with what we already have built:
http://desktop.arcgis.com/en/arcmap/10.3/manage-data/geodatabases/version-scenarios.htm#GUID-E09C2189-85FC-42AD-8DCA-2CC28B621CAF


texts: Newline appears after when converting shapefile to spatialite

texts:
I converted a layer from a shapefile format to spatialite. A line is added in the spatialite layer but record is not created in the attribute table (similar to multipart geometry). Shapefile to spatialite conversion was done just by right clicking the shapefile layer -> Save As spatialite.

The snapshot of the shapefile as well as spatialite layer is attached below.

I am curious to know what would be the reason?


texts: What is the difference between sentinelsat and sentinelhub?


texts: How to find out farthest north south east and west coordinates of a polygon?

texts:
I have hundreds of employees spread around the country. I need them to fill a form and they need to fill a mandatory decimal coordinates field.

They often write wrong coordinates that end up in the sea or neighbour countries and that messes up my analysis.

I want to write an Excel formula that prevents them from entering GPS coordinates that are not in the area of their responsiblity. In order to do that, I need to know the farthest north, south, east and west coordinates of a polygon, for example their city, province or village in ArcMap.

Could you please tell me how or point me to right direction?


texts: Converting values from st_area() to square meters

texts:
I am trying to find the area of a table in square meters. Note my spatial reference is in feet.
My query is something like :

 update counties set t_as = st_area(the_geom)/1609.34^2;

So basically i am trying to convert a value like 507863856.33256 to square meters. I want the area to be values to be something like 4.2361524 instead of 507863856.33256.

update: The SRID is 2263


texts: Dealing with gap between lidar-derived tree crowns that acquired from different las tiles using ArcGIS Desktop

texts:
I used lidR package and applied simple watershed algorithm for 2 las tiles. Between the segmented crowns there is a gap around 1 meters.

How can I deal with this gap?

I tried to "integrate" tool of ArcMap but all my crowns were integrated and this caused another problem.


texts:
I would like to use the Google tiles on LeafLet (RoadMap, Satellite,...).
I've tried this https://github.com/shramov/leaflet-plugins and it works fine, but I'm worried about the possible legal issues of this solution.
The following text is an extract of LeafLet FAQ.

The problem with Google is that its Terms of Use forbid any means of tile access other than through the Google Maps API.
You can add the Google Maps API as a Leaflet layer with a plugin. But note that the map experience will not be perfect, because Leaflet will just act as a proxy to the Google Maps JS engine, so you won't get all the performance and usability benefits of using Leaflet when the Google layer is on.

I think that's right but it depends on the point of view and the differences in interpretation of the Google Maps Terms of Use. For example Google says:

(a) No Access to Maps API(s) except through the Service. You must not access or use the Maps API(s) or any Content through any technology or means other than those provided in the Service, or through other explicitly authorized means Google may designate. For example, you must not access map tiles or imagery through interfaces or channels (including undocumented Google interfaces) other than the Maps API(s).

Here, LeafLet and this plugin could be consider as other technology.
Also Google says:

(g) No Use of Content without a Google Map. You must not use or display the Content without a corresponding Google map, unless you are explicitly permitted to do so in the Maps APIs Documentation, or through written permission from Google. In any event, you must not use or display the Content on or in conjunction with a non-Google map. For example, you must not use geocodes obtained through the Service in conjunction with a non-Google map. As another example, you must not display Street View imagery alongside a non-Google map, but you may display Street View imagery without a corresponding Google map because the Maps APIs Documentation explicitly permits you to do so.

At this point it could be interpreted as if we're not using a Google Map, but in a different point of view we could say we're using a Google Map under a LeafLet Map so strictly we're using it.
What do you think? Could I have some troubles with Google if I use this plugin?


texts: Where is street level data stored in mbtiles files?

texts:
I have a .mbtiles file loaded into an SQLite browser, trying to find street level data (road names, A/M roads etc) and I cannot figure out where any of it is kept.

I'm assuming it's stored in the tile_blob column, but I can't find any information on accessing the data, just how to render the image. Right now I don't need to render the map, just to be able to query it.

Can anyone offer some advice?


texts: How can I retrieve lat and lon rasters of a Geotiff or NetCDF image?

texts:
I am working with Landsat 8 and MODIS images for obtain Chlorophyll concentration

I processed Landsat images with ERDAS and get the chlorophyll concentration, (the product was a GeoTIFF image), I want to compare the information obtained with MODIS images in MatLab, for this I had to convert the GeoTIFF to NetCDF format with ArcGIS, but this did not work, because I need the raster of Lat, other raster of Lon and one raster of Chlorophyll in the same file.

I don't know how I can obtain these rasters.

Can you help me, please.


texts:
Cesium polylines follow the curvature of the earth, between two points, by default. How do I make arcs, which go above the ground, like in this example?

http://armsglobe.chromeexperiments.com/


texts: Points within "normal bounds" of a linestring in Python

texts:
I am working in Python (using Shapely). I have a list of Point(s) (the small circles in the image) and a LineString (the solid line with an arrow). I want to select only the points which are within the "normal bounds" of the LineString (the dashed lines show the normal bounds of the line - they intersect at 90 degress at the ends of the line).

I am using the below code to get the projection.

shapely_linestring.interpolate(shapely_linestring.project(shapely_point)

but if a point lies outside the "normal bounds", I still get a projection on the line at the nearest end. It does not give me a projection on the extended line, which I could have used.

Is there a better alternative to achieve identify such points?


texts:
I was selecting features using IFeatureSelection interface like as below.

IFeatureLayer2 featLayer = layer as IFeatureLayer2;
IFeatureSelection featureSelection = featLayer as IFeatureSelection;
featureSelection.SelectFeatures(spatialFilter, ESRI.ArcGIS.Carto.esriSelectionResultEnum.esriSelectionResultNew, false);

I selected all the features in the map as I used the spatial filter to the full map extent. But what I did not come across there are various enum options like
esriSelectionResultEnum.esriSelectionResultAdd. I was thinking that this add narrows down the search more like it would filter out more records as like this.

 //IQueryFilter2 qFilter = new QueryFilterClass();
                    //qFilter.WhereClause = "WORKORDERID != null";
                    //featureSelection.SelectFeatures(qFilter, ESRI.ArcGIS.Carto.esriSelectionResultEnum.esriSelectionResultAdd, false);

This immediate above snippet did not work for me But I am not very sure and I did not find much document on what scenario these options work and it would be great if it can be explained with snippet.


texts: Geoprocessing service saved data to scratch gdb

texts:
I created and published a geoprocessing service from a toolbox. It works with a file geodatabase on my local drive. When run as a tool in ModelBuilder, it worked well and wrote data into the file geodatabase perfectly. However, when I published it as geoprocessing service and triggered it via REST service, it wrote data into the scratch.gdb in the job directory on the ArcGIS Server directories.

How can I configure the tool/model, or service or arcgis server to make it work with the original file geodatabase instead of the scratch file geodatabase?


texts: Remove some islands from Europe json file

texts:
It's the first time I use a JSON file to create a map of Europe.
I downloaded these two files:

nuts0.json
nuts2.json

The first represents the states, while the second represents the regions of each state. I use them both because when you click on a country, the map is zoomed and I show the regions of the chosen state.

What I would like to do now is remove some islands that do not interest me (those circled in red).

How can I do that? Of course I want to remove them from both nuts0.json and nuts2.json.

UPDATE

These are the regions I would like to delete:

I was thinking of using mapshaper because it's online and it seems simple, but I can use any software to let me do what I want.

I didn't try to do anything because I don't know how to do, the only idea that comes to me is to edit the file by hand but this thing does not seem particularly intelligent, because I don't know how to change the arcs.

I would like to delete those circled. I have read somewhere that it is possible to establish a threshold, and to eliminate the smaller regions of this threshold. But there are regions (see Sardinia) that I care to keep and Sardinia is certainly smaller than the French Guaiana.


texts: Leaflet: Map focus issue after enabling a custom ZoomToArea control in Firefox

texts:
I'm trying to add a Zoom-to-Area control to Leaflet by extending the built in Control.Zoom plugin, but I'm fighting some "map-loses-focus" issue that affects only Firefox.

The function works as I intend in both Google Chrome and IE--once you click on the tool to enable, on mousedown you start drawing a rectangle, then on mouseup the action completes and the map zooms to the bounds of whatever rectangle you've drawn.

However in Firefox, after clicking/enabling the control, on mousedown it's like the map didn't have focus yet, and the dragging action (growing the rectangle) basically grabs a tile and just drags it around in the browser ..then when you mouseup the result is all "janked", it releases the dragging tile, still zooms to some portion of your rectangle, and fails to cleanup the rectangle in the map. BUT ..if you immediately do another zoom, it works as I would expect.

I have tried to programmatically force focus onto the map using two tricks, including one described here, but they don't solve the problem:

    // Prevent panning while drawing the zoom-to-area box.
    this._map.dragging.disable();

    // Try to return focus back to the map, necessary for FireFox.
    // 1..
    if (this._map) {
        this._map.getContainer().focus();
    }
    // 2..
    document.getElementById('map-canvas').focus();

I created a fiddle (same link as above) so you can see the complete code and compare the behavior in different browsers. To spot it quickly, I put my code extending the Zoom plugin between two comments:

    // start added
    //
       .. so my stuff would be here

    // end added
    //

I believe my _zoomToArea: function is where a work-around should likely appear.

Does anyone have a trick to prevent this unwanted behavior in FireFox?


texts:
I have a tempreture raster file, its min value is -125 and it's max is 50 I want to make a color-relief file of this raster,
I tried to use gdaldem to make it. It works but the problem is with its bounds, I used a color file text as bellow

-126 0 0 50
-83 10 10 100
-41 30 30 150
10 70 70 200
18.5 100 100 225
19 170 170 255
20 57 151 107
20.5 117 194 103
21 230 230 128
21.5 202 157 75
22 185 154 100
25 220 220 220
27 250 200 255
30 100 100 200

Well I think gdaldem first tries to make the image histogram equlizer and then make a color-relief, this is why my raster become as following

but it must be as

I tried to add -125 as nodata but nothing changed, what other options do I have to make a color-relief file from this raster?


texts:
Error "Dump Failed" minivolcado written in a temporary folder causes that my complement is closed which I am developing without finishing executing and closes QGis v.2.8.

What are these types of errors and what is happening so that it is generated?

Now uninstall QGis and delete all the folders and their registry, also delete the temporary files of the machine and reinstalled, but still giving the error.


texts: Seeking Python library for Sentinel SAFE files?

texts:
Is there a library for navigating Sentinel's SAFE product and tile archives (i. e. a nice object-oriented API for reading metadata, getting tiles, reading imagery files)

As far as I understand, sentinelsat only has utilities for searching and downloading the data. Sentinelhub has some code related to the structure of the SAFE format, but it is tailored to build SAFE archives from the AWS files and not general operations with SAFE.

Reading SAFE archives for analysis seems to be a common task, so just wanted to check that no available library exists before writing my own set of classes.


texts: Features are displayed as inline-variable

texts:
I run ArcGIS Pro 2.0 on windows 10. '

I run a model where the last tool is 'Sort' and the output is a featureclass to a geodatabase.

In my models I use the iterator inline variable, %Name%, from feature class iterator.

It names the files, as seen in windows file explorer, correctly, like: SortedPollution- where Pollution is %Name%

In map view the files are displayed as the output of the model, like: Sorted%Name%.shp.

Odd error, I am sure it has a simple solution.

In the case on the attached image you can see it displays the properties of "Sorted%Name%" when it is pointing to the file named correctly as "Sorteddkjord_v2".


texts: Showing map in Grass Layer Manager/Map Display via Python

texts:
I've created some raster(tif) maps using Grass libraries in Python. I want to opem them automatically in Grass GIS Map Display or Layer Manager after creation. I don't want to use GUI for this. What commands I need to use for this?

I'm using GRASS GIS 6.5.svn over Windows 7.

Thansk in advance..


texts: Does per-pixel viewing zenith angle (VZA) is identical between landsat sensors for a common path/row?

texts:
Due to the distance from the nadir viewing vector of the satellite, variation in viewing zenith angle (VZA) is observed from the center to the outer edges of the scene (as shown in figure one).

Because Landsat 5 TM, 7 ETM+ and 8 OLI all have a field of view of 15° (±7.5° from nadir) I was wondering if I can expect changes in pixel-based VZA over the entire Landsat time series for a given path/row, or if it is constant. And if it's constant, it is within sensors only or also between sensors ?

Furthermore, Is it expected to have a difference between bands for a same scene ?


texts:
I have the following code which works fine:
select vissoort, count(1), ST_Buffer(ST_GeomFromText('POINT(5.341248 51.615590)',4326):: geography, 2500)
from   visvangsten
where  st_intersects(visvangsten.locatie,
ST_Buffer(ST_GeomFromText('POINT(5.3412480 51.615590)',4326):: geography, 2500))
group by vissoort
order by 2 desc

Now I want the same function but then selecting the features within a polygon instead of the circle/buffer.
I tried things like this but nothing worked:
select vissoort, count(1), ST_asText( ST_Polygon('LINESTRING(5.303 51.629, 5.387 51.626, 5.393 51.588, 5.281 51.592)'::geometry, 4326) )
from   visvangsten
where  st_contains(ST_asText( ST_Polygon('LINESTRING(5.303 51.629, 5.387 51.626, 5.393 51.588, 5.281 51.592)'::geometry, 4326) ), visvangsten.locatie);
group by vissoort
order by 2 desc limit 1

The database table looks like this:
  id ([PK]bigint) | datum(date) | vissoort(character varying) | locatie(geometry)
------------------|-------------|-----------------------------|------------------
      15729       |  2007-06-23 |           Blankvoorn        |  0101000...etc.
------------------|-------------|-----------------------------|------------------
       etc.       |    etc.     |              etc.           |       etc.

Does someone know the answer?


texts:
I get this error when I'm trying to import a raster to RSAGA:

> rsaga.import.gdal(paste(rast_name, ".tif", sep=""),"RASTER", env= my_env)

Module 'GDAL: Import Raster' not found in SAGA library 'io_gdal'.
Check if module name has changed (or is misspelled)?
The following (non-interactive) modules currently exist in this SAGA library:

$io_gdal
NULL

However that file is found from modules, defined to my_env:

> my_env <- rsaga.env(workspace=substr(path_saga, 1, nchar(path_saga)-1),
+ path=sagapath, modules=sagamods)

> my_env
$workspace
[1] "C:/HY-Data/HRIIHIMA/data/DEM10_Fenno/"

$cmd
[1] "saga_cmd.exe"

$path
[1] "C:/Program Files/QGIS Wien/apps/saga/"

$modules
[1] "C:/Program Files/QGIS Wien/apps/saga/modules/"

$version
[1] "2.1.2"

$cores
[1] NA

$parallel
[1] FALSE

$lib.prefix
[1] ""

Screenshot:

Shouldn't this command work even without explicitly tell where the libraries area, like this:

> rsaga.get.libraries()
Error in dir(path, paste("^.*\\", dll, "$", sep = "")) :
invalid 'path' argument

It appears that rsaga can't find these libraries...

I'm using v. 3.1.2 R and v. 2.1.2 SAGA GIS

Any idea how to fix this?


texts: How to receive layer data features in geoserver?

texts:
I'm developing module in GeoServer that needs Layer feature members.
I can get it easily with WFS request like this

ows?service=WFS&version=1.0.0&request=GetFeature&typeName=workspace:layer&maxFeatures=50

But I need to achieve it internally in GeoServer(I mean something like this catalog.getFeatureMembersByName("workspace:layer"))

I know that I can send request in GeoServer internally but I don't think that is proper way. I don't want to talk directly with database because I need in this module GeoFence features.

Is it possible to achieve this?


texts: OpenLayers, WMS MinScaleDenominator and MaxScaleDenominator

texts:
In my OpenLayers map application I display many WMS layers from different sources. Some of them have display scale restrictions set by MinScaleDenominator and MaxScaleDenominator tags in WMS GetCapabilities.

What I want to achieve is to notify my application users to zoom in or zoom out to view particular layer.

How can I compare current map view zoom level to MinScaleDenominator and MaxScaleDenominator in OpenLayers?


texts:
I'm developing a PostGIS function to return a value when longitude and latitude are passed as variables to the function.

CREATE  FUNCTION flood_location(long float, lat float) RETURNS varchar(40) AS $BODY$
        BEGIN
                RETURN l.committee_new
                from land_ownership.council_owned_land l
                (long,lat),4236),l.wkb_geometry);
                where ST_WITHIN(ST_SetSRID(ST_MakePoint(lat,long),4236),l.wkb_geometry);
        END;
$BODY$ LANGUAGE plpgsql;

However the function is currently returning the following error when I trying to call the function.

select flood_location(-2.142,53.41862);

ERROR:  contains: Operation on mixed SRID geometries (MultiPolygon, 4326) != (Point, 4236)
    CONTEXT:  SQL statement "SELECT l.committee_new
                    from land_ownership.council_owned_land l
                    (lat,long),4236),l.wkb_geometry);
                    where ST_WITHIN(ST_SetSRID(ST_MakePoint(long,lat),4236),l.wkb_geometry)"
    PL/pgSQL function flood_location(double precision,double precision) line 3 at RETURN
    SQL state: XX000

The SRID of the land_ownership table is 4236?


texts: What is the difference between those two options of symbology export from QGIS?

texts:
What is the difference between those two variants of symbology export. I tried to export a point layer from QGIS to .tab format using both of them. "MapInfo Cartographic" True Type Font .ttf was used for symbology. The same font is used in MapInfo too.

Export settings:

But none of those two variants are successful.


texts: Why is my polyline on the wrong bearing?

texts:
I am trying to overlay duplicate segments of polylines on a map via javascript. I am using latlon.js and geo.js from here

The idea is that a user can input the distance they have travelled and when the map is refreshed, a second multipolyline will be added. To calculate the the point which represents where the user currently is, I am using a known previous location and known 'next location' to calculate a bearing between the two and then using that bearing, combined with the distance the user has travelled to calculate where they are now, before finally creating a new polyline based on that data.

The problem is that this new polyline is drawn incorrectly on my map and I don't understand why. You can see the problem at [removed]. Log in as User2/pass2 and then select 'challenges' from the nav bar. The most relevant code from map.js follows. I am using leaflet maps.

  var userDistance = GymloopChallenge.user_distance;
  var lastPoint;
  var comparisonArr = []; //used to calculate bearing
  if (userDistance <= dist1) {
      comparisonArr = [rioLL,laLL];
      distToAdd = Math.abs(userDistance);
      userPoly = [];
      lastPoint = rioLL;
      console.log('if 1');

  } else if(userDistance >= dist1 && userDistance <= dist2 ) {
      comparisonArr = [laLL,nyLL];
      distToAdd = Math.abs(dist2 - userDistance);
      userPoly = [line5];
      lastPoint = laLL;
      console.log('if 2');
  } else if (userDistance > dist2 && userDistance <= dist3) {
      comparisonArr = [nyLL,lonLL];
      distToAdd = Math.abs(userDistance - dist3);
      userPoly = [line5, line1];
      lastPoint = nyLL;
      console.log('if 3');
  } else if (userDistance > dist3 && userDistance <= dist4) {
      comparisonArr = [lonLL,sydLL];
      distToAdd = Math.abs(userDistance - dist4);
      userPoly = [line5,line1,line2];
      lastPoint = lonLL;
      console.log('if 4');
  } else if (userDistance > dist4 && userDistance <= dist5) {
      comparisonArr = [sydLL,tokLL];
      distToAdd = Math.abs(userDistance - dist5);
      userPoly = [line5,line1,line2,line3];
      lastPoint = sydLL;
      console.log('if 5');
  }

console.log("Distance to add: "+ distToAdd)
//console.log("COmparison array: "+comparisonArr);
  var bearing = comparisonArr[0].finalBearingTo(comparisonArr[1]);
  var userLocation = lastPoint.destinationPoint(bearing, distToAdd);
  //console.log(userLocation);

  userLineStartPoint = new L.LatLng(lastPoint._lat,lastPoint._lon );
  userLocationMapPoint = new L.LatLng(userLocation._lat,userLocation._lon );
  userPoly.push(new Array(userLineStartPoint ,userLocationMapPoint));
  console.log(userPoly);
    var userPolyline = new L.MultiPolyline(userPoly, {
      color: 'red',
      id:'poly2',
      weight: 3,
      opacity: 0.8

  });
  map.addLayer(userPolyline);


texts: What expression works with setLayoutProperty in mapbox

texts:
I am able to toggle the opacity of individual groups within a tileset. Unfortunately this does not hide the data when clicked on. As the layers are stacked this is a problem.
const toggleGroups = (layerName) => {
 if (mapRef.current && !loading) {
   let map = mapRef.current.getMap()
   let expression = ['match', ['get', `${layerName}_class`]]

   Object.keys(mapboxGroups[layerName]).forEach(key => {
     expression.push(parseInt(key), mapboxGroups[layerName][key] ? 0.58 : 0)
   })

   expression.push(0)
   map.setPaintProperty(layerName, 'fill-opacity', expression)
 }
}

The setLayoutProperty visibility works for other layers without groups and it also hides the layer data from being shown when layers are stacked.
I have tried to convert my function to work with setLayoutProperty but I get an error.
  const toggleGroups = (layerName) => {
if (mapRef.current && !loading) {
  let map = mapRef.current.getMap()
  let expression = ['match', ['get', `${layerName}_class`]]

  Object.keys(mapboxGroups[layerName]).forEach(key => {
    expression.push(parseInt(key), mapboxGroups[layerName][key] ? 'visible' : 'none')
  })

  map.setLayoutProperty(layerName, 'visibility', expression)
}

}
Error message reads : layers.risk.layout.visibility: expected one of [visible, none]
Is there another way of hiding the groups within a layer and also hiding their properties so when stacked the data is not visible when clicked?


texts: Interpolating Temperature over Elevation on ArcMap

texts: I need to overlay or create a new image file in ArcGIS that overlays temperature point data I have for California onto a 30m resolution DEM for that state. I have explored Co-Kriging, however I have not been able to firmly create a final product in which the temperature data is even remotely close to conformed to elevation. Any suggestions?


texts: Updating attribute table from an Excel file using Python

texts: I am kinda new to Python, please excuse for such a trivial questions. There is a shapefile (let's call it Materials.shp). I have an attribute table with different fields (for example, windows, doors, etc..). The script I am trying to write needs to update the values for those fields. The new values for those fields I get are from an Excel file. Is it possible to do that? What is the easiest way for the script to update the values of these cells? I am using ArcGIS for Desktop.


texts: Get ID from closest feature in different layer

texts:
ArcGIS Collector (Android)
ArcGIS Online

Scenario:
Users create and update features in a point FC via Collector.
Upon creation/update, I want Collector to:

Automatically get the ASSET_ID value from the closest feature in a line FC.
Insert the value from the line FC into an ASSET_ID field in the points FC.

What would it take to configure ArcGIS Collector to automatically perform a spatial operation like this?


texts:
I am trying to convert cursors written in 10.0 to use the new arcpy.da module in 10.2.1. If I step through the code in debug mode, it steps into the for loop; however, when I attempt to run the code I get an exception : iteration not started

code for 10.0

for row in updateCursor:
    if row.getValue("STREET_NUMBER_PREFIX") in ['N', 'S', 'W']:
        row.setValue("STREET_DIR_PREFIX",  row.getValue("STREET_NUMBER_PREFIX"))
        row.setValue("STREET_NUMBER_PREFIX", None)

code for 10.2.1

with arcpy.da.UpdateCursor(civicAddress, flds) as rows:
    for row in rows:
        if row[0] in ['N', 'S', 'W']:
            row[1] = row[0]
            row[0] = None

What can I do to solve my problem?

I have quite a few code blocks as above that deal with specific logic cases.

Here is the full code for that section, maybe it will help:

flds = ("STREET_NUMBER_PREFIX", "STREET_DIR_PREFIX", "STREET_NUMBER_SUFFIX", "CIVIC_ID")
flds2 = ("CIVIC_ID", "UNIT_NUMBER")
log.information ("Creating cursors...")
insertCursor =  arcpy.da.InsertCursor(subAddress, flds2)
#row = []

with arcpy.da.UpdateCursor(civicAddress, flds) as rows:
    for row in rows:
        if row[0] in ['N', 'S', 'W']:
            row[1] = row[0]
            row[0] = None

        if row[2] in ['N', 'S', 'W']:
            row[1] = row[2]
            row[2] = None

        if row[2] in unitLst:
            newrow = insertCursor.newRow()
            newrow[0] = row[3]
            newrow[1] = row[2]
            row[2] = None
            insertCursor.insertRow(newrow)

        if row[0] in unitLst:
            newrow = insertCursor.newRow()
            newrow[0] = row[3]
            newrow[1] = row[0]
            row[0] = None
            insertCursor.insertRow(newrow)

        try:
            if re.match("[0-9][A-Z]", row[0]) or re.match("[0-9][0-9][A-Z]", row[0]) or \
            re.match("[A-Z][0-9]", row[0]) or re.match("[A-Z][0-9][0-9]", row[0]):
                newrow = insertCursor.newRow()
                newrow[0] = row[3]
                newrow[1] = row[0]
                row.setValue("STREET_NUMBER_PREFIX", None)
                row[0] = None
                insertCursor.insertRow(newrow)
        except:
            pass

        try:
            if re.match("[0-9][A-Z]", row[2]) or re.match("[0-9][0-9][A-Z]", row[2]) or \
            re.match("[A-Z][0-9]", row[2]) or re.match("[A-Z][0-9][0-9]", row[2]):
                newrow = insertCursor.newRow()
                newrow[0] = row[3]
                newrow[1] = row[2]
                row[2] = None
                insertCursor.insertRow(newrow)
        except:
            pass

        try:
            if row[0].isdigit():
                newrow = insertCursor.newRow()
                newrow[0] = row[3]
                newrow[1] =  row[0]
                row.setValue("STREET_NUMBER_PREFIX", None)
                row[0] = None
                insertCursor.insertRow(newrow)
        except:
            pass

        try:
            #if row.getValue("STREET_NUMBER_SUFFIX").isdigit():
            if row[2].isdigit():
                newrow = insertCursor.newRow()
                #newrow.setValue("CIVIC_ID", row.getValue("CIVIC_ID"))
                newrow[0] = row[3]
                # newrow.setValue("UNIT_NUMBER", row.getValue("STREET_NUMBER_SUFFIX"))
                newrow[1] = row[2]
                #row.setValue("STREET_NUMBER_SUFFIX", None)
                row[2] = None
                insertCursor.insertRow(newrow)
        except:
            pass

    rows.updateRow(row)
    del row
    del newrow
    del rows


texts: Rpy2 Unable to Load in Fresh QGIS Install

texts:
after installing qgis the a message given below is shown

Unable to load the plugin: Unable to load required package rpy2.
Please ensure that both R, and the corresponding version of Rpy2 are correctly installed.

please let me know why this is fault


texts: How to calculate X and Y coordinate of a surfacic layer in QGIS

texts:
I have a surfacic layer in QGIS 2.12.3. I need to calculate the X and Y coordinates for the centroid of my polygons. I created a new field in my attribute table. When I try to use the field calculator, it does not work. In order to solve the problem, I created a new layer with the centroïds of my polygons. In this layer I was able to dertermine the X and Y coordinates. Now, I want to transfer this information to my orignial surfacic layer.

Is it possible to calculate the X and Y coordinates directly from a surfacic layer? If it is possible, how could I do it ? If it is not possible, how can I transfer the X and Y coordinates from my centroïd layer to my surfacic layer attribute table?

I used this post to calculate the X and Y coordinates for my centroïds: How to write a expression in field calculator for finding $X and $Y?


texts: Heuristic relationship between number of satellites and GPS precision?

texts:
I'm a statistician conducting an analysis on animal movement based on data collected using ATS Wildlink 500 collars. As part of this analysis, I'm estimating GPS error assuming that the error in the measurement is normally distributed with the mean being the true location (x and y) with some standard deviation that is equal in the x and y dimension. I'm ignoring accuracy at this point, or rather, I'm assuming that bias is zero.

Obviously, there is way more to this, but as a first cut analysis, I'm happy with this assumption. However, the one thing that I want to capture at this stage is the likely degradation in precision based on the number of satellites associated with a given GPS fix. What I'm looking for is a simple relationship between the number of satellites and the standard deviation of the GPS error given that the fix was successful. Less than a minimum number of satellites (4) and the error is infinite (i.e. the fix fails). So really I'm just looking for the relationship when the number of satellites is above this number (say between 4 and 12).

Is it linear?
Sigma = A + B * X, where is  X is the number of satellites?

Exponential? Something like:
Sigma = exp(A + B*X)

Unfortunately, I don't have raw data where I could investigate this relationship myself and I'm having trouble finding good resources online. So I'm going to have to fit this relationship blind to a certain extent. I'm just looking for a general guideline on how precision improves with increasing number of satellites.


texts:
I am trying to mosaic 3 S2 L1C images. I used the multi-size mosaic to prevent resampling as preprocess for the mosaic-only. However the resulting image has 0 values.
This is the result:

Since S2 images were big, I also tried subsetting the 3 images and repeated the steps above. The result is just the same. How to do this correctly? Should I just use resample then mosaic (not multi-size mosaic)?
Note: I also posted this as question in SNAP forum but it is still for approval.


texts: Creating polygon with extents of raster file in QGIS

texts:
I have a tiff image file and want to get a polygon with the extent of the tiff.

Is there a direct way of doing this?

Right now I have to convert the image to vector (with gdal_polygonize) and then dissolve the polygons.


texts: How to fix SRTM half-pixel shift with GeoServer

texts:
The SRTM data has an half-pixel shift (I assume that’s what is described in the last FAQ, and that GeoServer uses pixel corner origin ?  I couldn’t find this information).
Despite this shift being apparently widely known, I couldn’t find directions on the best way to correct and publish these data in GeoServer.
I can think of two options :

reprocessing the data files before importing them ;
defining a custom coordinate operation in GeoServer.

What is the recommended option ?  I guess the answer depends on the frequency of the requests to the server (the former option being preferred for high frequency requests to lighten the server computation needs), but maybe I am missing a third option ?
Data are loaded into GeoServer through ImageMosaic plugin, and here is what gdal_info has to say about the .tif files :
$ gdalinfo srtm_38_03.tif
Driver: GTiff/GeoTIFF
Files: srtm_38_03.tif
       srtm_38_03.tfw
Size is 6000, 6000
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (5.000000000000000,50.000000000000000)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   5.0000000,  50.0000000) (  5d 0' 0.00"E, 50d 0' 0.00"N)
Lower Left  (   5.0000000,  45.0000000) (  5d 0' 0.00"E, 45d 0' 0.00"N)
Upper Right (  10.0000000,  50.0000000) ( 10d 0' 0.00"E, 50d 0' 0.00"N)
Lower Right (  10.0000000,  45.0000000) ( 10d 0' 0.00"E, 45d 0' 0.00"N)
Center      (   7.5000000,  47.5000000) (  7d30' 0.00"E, 47d30' 0.00"N)
Band 1 Block=6000x1 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

And the .twf file :
$ cat srtm_38_03.tfw
            0.0008333333
            0.0000000000
            0.0000000000
           -0.0008333333
            5.0004166667
           49.9995833333


texts:
How does one access data from an Esri file geodatabase in Quantum GIS?

I've seen snippets of conversation that indicate it's possible with qgis 1.7+ and gdal 1.9+ but no clear recipe.

I recognize that in a few months or a year this will likely be a non question ("just click [add data] button"), but today, I'm stumped. ;-)


texts: Customizing editor widget for feature layer

texts:
Here is the code for the editor widget:

<esri:EditorWidget x:Name="MyEditorWidget" Map="{Binding ElementName=MyMap}"
                                      Width="300"
                                      AutoSelect="False"
                                      GeometryServiceUrl="http://tasks.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer"
                                      ShowAttributesOnAdd="True"  />

This is provided by Esri. Some of these buttons have no use for our clients.

Is there any possible way to get rid of them?

Also, can I use some of the features from the editor without actually having the widget, such as the Add?

I have not seen any documentation for this.

Here is a picture of what the editor widget looks like:

The only buttons I am interested in are the Save Edits, Add, and Clear Selection.


texts: No Postgis after installing opengeo suite

texts:
After the installation of the opengeo suite; everything seems to work fine, except for the postgis extension of postgresql. In the dashboard all the postgis options are greyed out.

The system is running Ubuntu 12.04 (Precise) LTS Bare OS.

Is there a way to install postgis afterwards with the templates from opengeo?


texts: QGIS export selection to shapefiles is extremely slow; using only 5 to 10% CPU

texts: I am using QGIS 3.2.3. I am trying to do some simple actions. Selecting by attribute and exporting to shapefiles. I am using a OSM file as a source which is barely 60MB. But for some reason the process is incredibly slow (as in, it runs for twenty minutes and only reaches 5% of the export). When I look at my task manager QGIS is only using between 5-10% I suspect this is not much, but I don't know how to increase it (or even if I should.) Should I make it use more CPU? Or should I find another way to speed it up?


texts: Light normalization of Orthophotos to increase comparability

texts:
I would like to work with the RGB-values of different orthophotos. Now I noticed that the light conditions are not the same in the different orthophotos. I was wondering if there is a sound way to make these images a bit more comparable. I found quite some studies about 'multitemporal normalization' for exactly the same region. My images are not images of the exact same region though, but they are close by each other.

If I adapt the min-max values I might falsify too much because different objects in the photos reflect differently of course.

Does anyone have a sound solution for adapting these images (at least a bit). Or do you have recommendations what I could search for (because I run out of search terms)?


texts: GeoAlchemy get closest point to another point

texts:
I have loaded into my db the cities1000.txt available here. I am attempting to create a query using GeoAlchemy to return the nearest city when I pass in Lat and Long. I currently have the query setup using ST_ClosestPoint, but it returns the same city no matter what lat/lon values I pass in. Here is the query I have written so far:

q = db.session.query(cities1000,
                     func.ST_ClosestPoint(cities1000.c.geom,
                                          func.Geometry(func.ST_GeographyFromText(
                                              'POINT({} {})'.format(lon, lat))))).first()

Is ST_ClosestPoint what i'm looking for? Performance is a concern.


texts: show polygon from database in browser in geodjango

texts:
I tried to retrieve the polygon which I saved into postgis database. (all polygons)
Here is my code:

models.py

from django.contrib.gis.db import models

class MyMap(models.Model):
    name = models.CharField(max_length=50, blank=True, null=True)
    timestamp = models.DateTimeField(auto_now_add=True, auto_now=False)
    geom = models.GeometryCollectionField(srid=4326, null=True)
    objects = models.GeoManager()

    def __unicode__(self):
        return self.name

forms.py

from django.contrib.gis import forms
from .models import MyMap

class MyMapForm(forms.ModelForm):

    class Meta:
        model = MyMap
        fields = ['geom']

views.py

from django.shortcuts import render
from django.http import HttpResponse
from .forms import MyMapForm
from .models import MyMap

def map(request):
    map_info = MyMap.objects.all()

    map_data = {
        "map_detail": map_info
    }

    return render(request, 'mymap.html', map_data)

def ajax(request):
    if request.POST.has_key('geom'):
        form = MyMapForm(request.POST or None)
        ...

mymap.html

<form action="" method="POST">
    {%csrf_token%}

 // I use google map api
</form>

    {% for details in map_detail %}
        {{detail.geom}}
    {% endfor %}

I don't know where is my problem which I can not show the polygons in the html file !


texts: QGIS 3.4.2 not loading .tif files correctly

texts:
I have a surface which has been exported from Civil3D and added to both QGIS 2.18 and 3.4.2

When added they both look very different and the imported image and values in 3.4.2 is not very helpful in that I cannot see what the null values of the image are to remove them.

This cannot be corrected by changing the max. and min. values for the layer either as it just changes the color tone for the whole of the image as if it is all just the one height throughout.

Has anyone else experienced this when adding .tif files to QGIS3.4.2 or is there a new setting somewhere that I have not switched on or set correctly?

View in QGIS 3.4.2:

View in QGIS 2.18.21:


texts: Options for manually identifying road defects from Terrestrial Lidar

texts:
I am trying to create a process to identify road defects manually from Lidar. We have Erdas 2013, ArcGIS10.1 and Global Mappper 14.2 as options for paid software.

So far Merrick Mars displays the data the best but (the free version) does not allow one to draw polygons.

With ESRI the best I could get is

Unfortunately you can't zoom any further to view defects.

With ERDAS the best I could get is

With Global Mapper the best I could get is

I can manually assign a color to each range of elevation but this is very cumbersome.

So any ideas on what I am doing wrong? how can I get a similar view to Mars in any of the other software?

I need to then create a polygon around each road defect.


texts: Interactively determining upslope and downslope areas on DEM?

texts:
I know SAGA GIS can be used to interactively generate maps of upslope and downslope areas for a cell within a DEM, but the output is pretty ugly.

An upslope area is defined as all cells (for a particular routing algorithm) that contribute runoff to the cell of interest.

A downslope area is defined as all cells (for a particular routing algorithm) that receive runoff from the cell of interest.

Are there any other tools that perform the same calculation but produce nicer results?

Here's a acreenshot of the DEM in SAGA GIS, followed by the output from using the  Upslope Area [interactive] tool.  As you can tell from the two images, it becomes difficult to keep track (visually, which is all I hope to do) of the general location of the generated upslope area within the original DEM image.

Efforts in SAGA GIS:

Using the Upslope Area [interactive] and Downslope Area [interactive] tools, with the D8 (SFD) algorithm, I was able to generate the results below:

Upslope area (also known as contributing area) is in yellow; downslope area (also known as dispersal area) is in green. The point of interest is framed in a red box. Generating this image required running each tool on the DEM, then turning any non-area cell (with a 0 value instead of a value of 100 for the area cells) into a NODATA cell using the Reclassify Grid Values tool. That allowed layer of all the results together.

Efforts in Whitebox GAT:

Using the Trace Downslope Flowpaths tool, I was able to trace downslope and upslope flowpaths using the D8 algorithm.  I created a D8 Flow Pointer for my regular DEM using the D8 Flow Pointer tool under Flow Pointers under Hydrological Tools.  Then I created a D8 Flow Pointer for an inverse DEM (DEM values multiplied by -1, used the Multiply tool under Mathematical Analysis).  Then I created a point as using the instructions @WhiteBoxDev provides in the comment below - creating a new shapefile and adding digitizing a point on the DEM as a feature in the shapefile.  I then than ran the Trace Downslope Flowpaths tool using the appropriate DEM (regular or inverse) and the newly created point file.

The same point was NOT used between the two programs.  The main difference in the results, the hydrological validity/accuracy of which is certainly up to debate, is that the upslope area in SAGA GIS shows how many points can feed into the same cell, while the WhiteBox GAT upslope area results show single line flowpaths from the inverse DEM.  This is due to the difference in the use of the two available tools.  SAGA's upslope tool determines which upslope cells could eventually feed into the downslope cell of interest, allowing for an "area" to be generated.  CIn Whitebox GAT, we calculate a "reversed" downslope flowpath from an inverse DEM, generating one (possibly the steepest?) upslope flowpath, which shows the path of cells that runoff from the cell of interest would flow through with the inversed elevation values.

Once again, the hydrological accuracy of any of these results is certainly up for debate.


texts: Exporting Satellite Imagery from ArcGIS Desktop for Global Mapper?

texts:
I have satelite imagery in ArcGIS and I want to open it in Global Mapper.
The problem is that the after exporting the image it looks like the color of the image is changing.

Here, below I have the image that I want to export from ArcGIS :

Here is the way that I export the file :
1.Right click on the layer and export data.

Here,below is the same file after i opened it in Global Mapper :

So, any idea how to fix it?


texts: Efficiently Searching MXDs for Folder Connection

texts:
We are currently trying to restructure some of the companies GIS data. It is located on an internal network drive (no SDEs). For each year, there are 1000+ MXDs with map layouts that still need to be working after the restucturing process. Therefore im trying to analyse our MXDs before we start changing the data links using the .findAndReplaceWorkspacePath() function.

What I've done so far: 1. For each year, I've listed all MXDs available (very fast, using os.walk) and 2. In a second loop, I now want to find all MXDs that contain data pointing to a certain directory. The search process involves arcpy.mapping, which I find terribly slow... Here is my function for the search process:

def MXDchecker(inList, searchString):
    internalList = []
    internalCount = 0
    counter = 0
    for mxdFile in inList:
        print(u"\tChecking MXD: {}".format(mxdFile))
        try:
            mxd = m.MapDocument(mxdFile)
            # create list for links to v-drive:
            shortList = []
            # start search, find all links to v-drive
            for df in m.ListDataFrames(mxd):
                for lyr in m.ListLayers(mxd, data_frame=df):
                    if lyr.supports("DATASOURCE"):
                        if searchString in lyr.dataSource:
                            shortList.append(lyr.dataSource)
                            print(u"\t\tFound v-link.")
            # Check if list has entries:
            if shortList:
                internalList.append([mxdFile, [shortList]])
                internalCount += 1
            # delete mxd
            del mxd
        except:
            print(u"\tError for processing MXD {}".format(mxdFile))
        if internalCount == 100:
            counter += internalCount
            internalCount = 0
            print(u"\tFound {} MXDs with v-links...".format(unicode(counter)))
    counter += internalCount
    print(u"Found a total of {} MXDs with v-links.".format(unicode(counter)))
    return internalList

The returned list contains all relevant connections and the mxd path. I tested this on only 5 MXDs, some of which contain multiple dataframes and numerous layers. The processing times are very long (> 10 min), for a run on 1000 + MXDs it would take a full day at least, without even changing any of the data links...

My questions:

How can I speed up this process in general? Is it better to first find all MXDs containing any connection to the v-drive and then iterate on them, straight with the processing?

Note that:

I cannot copy the MXDs to a local drive for processing, it has to be done in the given location
I tried using the new ArcGIS Pro arcpy.mp module, but I found the structure to be too different for the use with "old" MXD files.
I'm working with Python 2.7 for ArcGIS 10.6.1 with an Advanced License. I do also have access to ArcGIS Pro (Advanced).


texts: From pgRouting shortest_path back to osm node id

texts:
I'm successfully performing path finding queries like shortest_path on data imported with osm2pgrouting.

The tables I have inside my pgrouting DB are the following:

               List of relations
 Schema |        Name         |   Type   | Owner
--------+---------------------+----------+-------
 public | classes             | table    | user
 public | geography_columns   | view     | user
 public | geometry_columns    | table    | user
 public | nodes               | table    | user
 public | spatial_ref_sys     | table    | user
 public | types               | table    | user
 public | vertices_tmp        | table    | user
 public | vertices_tmp_id_seq | sequence | user
 public | ways                | table    | user

I get the way of a given path between a source and dest node this way:

SELECT * FROM ways
   JOIN
    (SELECT * FROM
       shortest_path('SELECT gid AS id, source::int4 AS source, target::int4 AS target, length::double precision AS cost FROM ways',3,7,false,false))
   AS route ON ways.gid=route.edge_id;

ways's table has the following structure:

 Table "public.ways"
    Column    |       Type       | Modifiers
--------------+------------------+-----------
 gid          | integer          | not null
 class_id     | integer          |
 length       | double precision |
 name         | character(200)   |
 x1           | double precision |
 y1           | double precision |
 x2           | double precision |
 y2           | double precision |
 reverse_cost | double precision |
 rule         | text             |
 to_cost      | double precision |
 the_geom     | geometry         |
 source       | integer          |
 target       | integer          |

Now gid, source and target values aren't original osm node's ids (the node's ids contained in the .osm file from which the table has been filled). I suppose that pgRouting create custom ids in order to index ways and edges.
From what I can understand the original ids of nodes contained in the .osm files are the ones listed into nodes table.

I would like, given a row of ways table, recover the original osm id of source and target.

How can I do that?is it possible?


texts: MapBasic - Create a Table array based on a variable

texts:
I wanted to create a temporary table as per my variable.

If i have a value of 10, it will create a table with 10 rows.

How should i start my code in MapBasic?


texts: Styles don't get Tag when imported into QGIS

texts:
I have a .XML file with lines and fills as well as .svg Images that I want to use as markers. When I Import them into QGIS I want all of the styles to get the tag "IBF". But this doesn't work and I have no idea why. I can Import the styles, they are under the tag "all" and I can manually add each element the IBF tag and then the style appears also under the IBF tag, but why isn't it automatically under the IBF tag?. Is something wrong with my Settings? Here is a Picture of my Settings.

I'm in the "IBF" Tag when I click on Import and I wrote "IBF" in the "additional Tags" field.

Any ideas on what I missed?


texts: Exporting QGIS layer as shape through REST API

texts:
I would like to export a layer described in QGIS Server project file as Shape through REST API. For example GeoServer's GetFeature accepts outputFormat parameter value shape-zip which in turn returns Shape zip archive for the desired GeoServer layer.

Is there a way to accomplish the same with QGIS Server? GetCapabilities on my QGIS instance declared that GetFeature has following output formats:

<ows:Parameter name="outputFormat">
    <ows:Value>text/xml; subtype=gml/2.1.2</ows:Value>
    <ows:Value>text/xml; subtype=gml/3.1.1</ows:Value>
    <ows:Value>application/vnd.geo+json</ows:Value>
</ows:Parameter>

So I was wondering maybe I have missed something and maybe there is a way export layer data in same manner as with GeoServer (not looking for a plugin)?


texts: Large dataset processing via cloud computing

texts:
I need to process a large dataset.  Aside from processing locally on my computer, can I install, for instance, arcgis10 on a cloud service like amazon.  I don't have esri server, only arcgis10.

I'm not interested in storage or distribution, just processing.

I have ample bandwidth so upload/download isn't an issue.

Might seem like a silly question, but thought it was worth asking.


texts: Attribute table: Add columns from one layer to the other one (both shp types)

texts:
I have two layers, both shp type. One contains villages with all their details, the other one is a result from a point sampling tool, so two columns - population density and elevation.

What I would like to do is to add the two columns to the attribute table of villages layer.

I went to the attribute table of villages layer and tried to add these two columns but without result.

Is it possible? If yes, how can I do that?


texts: How to view locally generated XYZ map tiles?

texts:
Is there an easy way to test (view on a map) my locally generated XYZ map tiles?

I've used OpenLayers and Leaflet to display non-local XYZ map tiles, and I could conceivably push these tiles to a remote server and test them that way.  It just seems to me there should be an easier way...


texts: Why do layers appear in the wrong order in QGIS 2.2?

texts:
I have a problem with in which order my vector and raster layers appear in QGIS. I have raster that I want to be the background of my other layers (vector). I place the vector layers on top in the TOC, but some of them still end up under the raster layer in the display?
I have check there projection and that's the same for all layers. I have also tried putting the raster layer on top in the TOC, but the resault is still the same.

I just want to show a few vector layers on top of a raster, but in this version something goes wrong. I have never had any trouble with this in earlier version of QGIS.

Any ideas on how to solve this?


texts: How to symbolize null values in ArcGIS using Quantities symbology?

texts:
One numeric field in my data involves null values (missing data). When I specify it as quantities in symbology, it turns out that the polygons associated with null values are taken out of the map. How can I represent these values (i.e., their polygons) on the map? I read some where to convert null values to 0, but doing so is misleading as these are missing entries, not 0.

I get into this situation after joining the shapefile with an attribute table, but I guess it happens with any data that includes missing data. Is that so, or is only unique to the joining procedure?


texts: Choosing LiDAR tiles with vector data using LAStools?

texts:
So the problem is that I have a huge amount of LiDAR tiles (in .laz format), but need only part of them that are within the borders of the city.

When I load all the tiles to LASTools, I have something like on the photo below:

The blue line is the border, everything outside doesn't need to be used.

Is there any LASTools tool that would help me choosing tiles I need?


texts:
I am running a process in FME 2014 to translate some data from an Excel to a PostgreSQL database.
My created table needs to have two primary keys (thus a composite key).
I have selected from the writer properties the columns which I want them to be PK but when I run the translation I get:

"Multiple primary keys specified in definition for table 'bom'. Only one primary key can be specified per table. Discarding primary key request for column 'item_code'"

How can I define composite keys in FME? Is it possible at all? If not, what are the alternatives?


texts:
I'm interested in understanding the relative merits of using GeoJSON or shapefile formats for handling geo-spatial data, particularly in an open-source, Python (CKAN project), web-development context.

GeoJSON is well supported in the CKAN data management python project (through Open Layers), but our clients seem to work with shapefiles.  I'm trying to understand whether I should just ask clients to run their shapefiles through a simple GeoJSON convertor (perhaps even implement this myself upon file upload), or whether there is a good reason to keep their data stored in shpaefile format.

Struggling to find a decent technical comparison of the file formats.  Anyone care to break it down for me?


texts: Error while exporting GeoJSON from PostGIS

texts:
$ ogr2ogr -f "GeoJSON" test.geojson -t_srs EPSG:4326 PG:"host=localhost user=postgres dbname=sp_join port=5432" -sql 'select gid, geom from test'

I am using the following code in ogr2ogr to export a table from PostGIS to GeoJSON, but get the following errors:

ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported
ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported
ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported
ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported
ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported
ERROR 1: Reading EWKB with 4-dimensional coordinates (XYZM) is not supported

How can I avoid 4-dimensional coordinates?


texts: Generate MaxScaleDenominator in geoserver

texts:
My geoserver is providing few overlayers with specific informations such as poi's or fire-hydrants. I'm using scale dependent styling, so it's not rendering this features above specified scale, but it still tries to render it and returns empty tiles. I've found WMS 1.3.0 includes MaxScaleDenominator and MinScaleDenominator tags, which can be provided in the getCapabilities result to show the optimum visualization ranges for each layer. I want to use it to learn my application not to ask about tiles above specified scale (cause they're empty there) so it will improve performance and reduce bandwidth usage.

I searched all geoserver properties (general, layer, caching) and I was not able to find any place to set this feature up, so I don't think there is some place in web-gui (or there is??), but maybe there is some resolution to write it somewhere in configuration files?

This is some link to my geoserver, e.g. layer is osm:Hydranty_active

After some investigation I've found that in layer osm:Hydranty_active served by WMS this information is included, but the same layer served via gwc don't have it, so trouble will be with GWC not Geoserver


texts: Circumventing QgsCredentials script called from Pre-execution?

texts:
I've got a stored project that, when opened is connecting to my PostgreSQL DB, and is loading several layers. All good.

I'm using a script that is called from Pre-execution (Processing/Options and configuration)
that is using the QgsCredentialDialog to distinguish between admin and regular user.

def LoginCredentials():
  uri = QgsDataSourceURI()
  uri.setConnection("localhost", "5432", "Blah", None, None)
  connInfo = uri.connectionInfo()

  (success, user, passwd) = QgsCredentialDialog.instance().get(connInfo, None, Nones)
  if success:
    uri.setPassword(passwd)
    uri.setUsername(user)
    iface.messageBar().pushMessage("Welcome!", "Welcome "+str(user), level=QgsMessageBar.INFO, duration=3)

def run_script(iface):
  LoginCredentials()

Issues:

one of my 7 layers is able to be loaded before the credentials dialog pops up.
the credentials dialog is printing a string just above the username/password textboxes with dbname, host, port, user, password and sslmode, all filled out! (I'm guessing the values are read from the .qgs project file, but I rather it didn't show it)
if I keep hitting the Cancel button I can see that the layers are loaded one by one until everything is loaded and the credentials dialog just disappear.

I'm very new to QGIS (4 days), so it might be that I've overlooked something very essential.


texts: Not able to pass the parameters to SQL view in GeoServer

texts:
I created a sqlview(villagefilter_view) in GeoServer and added one parameter %vilselected%

select * from niphad_survey where village=%vilselected%

Name: vilselected
Default value: 'Ranwad'
Validation regular expression: /^[a-zA-Z]+$/

I tried to pass the parameter from a selec option onchange value with vilSelectForView name.

Following is my code to genarate request url from GeoServer

var niphadSurveyVectorSource = new ol.source.Vector({
        loader: function(extent) {
            $.ajax('http://localhost:8080/geoserver/postgis_22_sample/wfs', {
                type: 'GET',
                data:{
                    service: 'WFS',
                    version: '1.1.0',
                    request: 'GetFeature',
                    typename: 'villagefilter_view',
                    viewparams: "vilselected:"+vilSelectForView,
                    srsname: 'EPSG:3857',
                    bbox: extent.join(',') + ',EPSG:3857'
                }
            }).done(function(response) {

         niphadSurveyVectorSource.addFeatures(formatWFS.readFeatures(response));
            });
        },
        strategy: ol.loadingstrategy.bbox,
        projection: 'EPSG:3857'
    });

After page lode the default parameter data is loading on map but onchange my select option the parameter is not passing to request URL.

How can I load the dynamic view data on map by onchange with the village name?


texts: Ordering polygon neighbors in clockwise order using ArcMap?

texts:
In ArcMap, I am using the Polygon Neighbors tool to generate a list of each polygon's adjacent polygons, or 'neighbors.' Currently, ArcMap orders the polygons using each polygon's numeric ID. What I would like to do is order each polygon's neighbor list in clockwise order relative to the center of the "source" polygon.

For example, in the picture below, if Polygon 1 is my source precinct, the neighbors currently appear ordered 2,3,4,5,6,7. I would like the ordering to be clockwise - say, 6,3,7,5,4,2. (The start of the list does not matter, only that they are clockwise relative to the starting neighbor in the list.)

Is there a way to do this in ArcMap?


texts: What is the most efficient way of displaying GPS in real time using OpenLayers?

texts:
i'm would like to display GPS traces (not only a position) in realtime on OpenLayers. What is the most efficient way to do it? Sending KML or GPX files to browser or sending data via GeoJSON?

br
d


texts: How to write a QGIS Action which opens a PDF from a server?

texts:
I need to open a PDF from an action in QGis.
The PDF rests on a server, not local.

So far I tried this (as Python):

from PyQt4.QtCore import QUrl;
from PyQt4.QtWebKit import QWebView;
myWV = QWebView(None);
myWV.load(QUrl(('https://myserver.com/[% "path" %]));
myWV.show()

This just opens an empty window. It works with a png.

Can anybody tell me how to do that.
Thanks!


texts: Focal operations / aggregations on stars_proxy object in R

texts:
I have many NetCDF-files that I read in as stars_proxy-object. I also have a geopackage that is a multipolygon. I Would not like to extract the mean value for each polygon for each stars_proxy-object.
I don't event know if a stars_proxy object is the best way to do so. But the stars-package provides such great functionality with st_apply to aggregate over time or spatially, that I just used the stars-package to read the data. There is the st_extract-function, however this is only for points.
I'm searching for something like the raster::extract or exact_extract-functions that also aggregate cells-values intersecting with polygons.
Maybe a "normal" raster is better. I'm just not really sure and would be happy for any pointer:)


texts: gdal2tiles_parallel only creating one file per zoom level

texts:
I'm trying to create a TMS layer from a gdal VRT dataset. To speed thing up, I'm trying to use gdal2tiles_parallel (https://github.com/GitHubRGI/geopackage-python/wiki/Usage-Instructions-for-gdal2tiles_parallel.py). The command I'm using is:

gdal2tiles_parallel -e -z '5-15' --format=PNG sp2014_b321.vrt sp2014_b321_tms/

The program starts fine, the directory is created and the Z (zoom) folders are also created. However, only one X folder is created inside each Z folder and only one png is created for each zoom level. The program ends without error but the TMS is incomplete. Has anyone seen this? Any idea how to fix is?


texts: Converting 13 band .tif images to RGB using QGIS?

texts:
I have sentinel 2 images. Each image has 13 bands (bands are not individual images). I wanted to convert them to only RGB images and store them in a file to use for another purpose.

I have tried

raster - Conversion - PCT to RGB

But it does not work. this is the error:

 File "/usr/bin/pct2rgb.py", line 121, in

 ct_size = ct.GetCount()

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

is there a fix to this or any other methods to achieve converting multiband to rgb.


texts: Overlaying rectangular polygon grid to extract tiled map subsets?

texts:
I want to be able to create/extract subsets(or tiles) of map data to new layers so that I can export them as tileable jpegs.  The tiling is important for my project. My strategy is to place a grid of rectangle Polygons on individual layers and use them to extract map subsets from an underlying Shapefile.

When I create a new project the CRS is set to WGS-84 EPSG-4326.

I have an underlying map - a Shapefile of building plots - when I import it, its CRS is set to CRS(EPSG:4326, WGS 84). My area is in Malmö, Sweden.

I then use the plugin quickWKT and paste a WKT file of a set of polygons: the file looks like this:

POLYGON((1449008 6180085,1449008 6163387,1414376 6163387,1414376 6180085,1449008 6180085))

etc etc...

(I have 9 polygons, 3x3).

The new layer created from these polygons is a shapefile layer (I think) and has same CRS as the project CRS and underlying layer. Now from what I understand, the coordinates used for these polygons are system RT90 and not WGS. The polygons do not map onto my underlying layer of house plots. They are far off in the distance.

Can anyone help me in a direction?


texts: How do I query in natural=strait from OSM databse?

texts:
I'm trying to query in straits, fjords, and bays from a PostGIS Openstreemap database. The straits I'm looking for are tagged natural=bay, natural=strait, and natural=fjord.

I used the following code:

(SELECT
        way, natural
      FROM planet_osm_line
      WHERE natural IN ('bay', 'strait', 'fjord')
    ) AS straits

But, it returns:

 Postgis Plugin: ERROR:  syntax error at or near "FROM"
LINE 3:           FROM planet_osm_line
                  ^

in executeQuery Full sql was: 'SELECT * FROM (SELECT
            way, natural
          FROM planet_osm_line
          WHERE natural IN ('bay', 'strait', 'fjord')
        ) AS straits LIMIT 0'

I am able to query docks / piers. And the SQL I used is very similar so I can't understand what I'm doing wrong.

This is what I used for docks and it works:

(SELECT
        way, man_made
      FROM planet_osm_line
      WHERE man_made IN ('pier')
    ) AS docks


texts: ArcMap ModelBuilder iterate Raster and Feature Class with Raster Calculator

texts:
In our model we calculate raster datasets with the cost distance and the raster calculator. We want to multiply these with a field value from the table testbezirk. We have already tried a few things, no error message comes up or ArcMap throws exception and there is no message, but all rasters are multiplied with the last value from the table testbezirk. It should be so that the field value (long variable) is calculated from the same row from which the cost distance is calculated. So :

The Test District table has the source items for the cost distance and then the field value Cyclist is to be taken from the same row and multiplied by the raster data p_%Name%.
We are not getting anywhere and have already tried the whole thing with sub models, but it does not work. Does anyone have an idea?
In detail: the main model contains the sub model cost distance and probability. In the cost distance the feature classes "testbezirk" are iterated with a source raster. The datasets are stored in a geodatabase named "costdistance" by name cd_%Name%. In the probability submodel the costdistance datasets are iterated from the "costdistance" geodatabase and calculated using the raster calculator. These raster datasets are then stored in another geodatabase "probability" by name p_cd_%Name%.
The datasets from the geodatabase "probability" should now be iterated and multiplied with the corresponding field value "cyclist" from the table "testbezirk", so the field value should also be iterated. Whenever a raster is iterated, the corresponding field value is to be iterated.
Example: r_p_cd_1 * 66, r_p_cd_2 *51, etc...


texts: GDAL NULL Path Value On Production Server

texts:
I am using c# and visual studio.

I have an app using GDAL that works perfectly on the Dev Server. When I transfer it to Production the app fails. I receive the following error message:

An error occurred while warping the file: System.ArgumentNullException: Value cannot be null.
Parameter name: path 1

   at System.IO.Path.Combine(String path1, String path2)

   at GDAL.Services.GDALService.Warp(String fileName, String destinationFileName, String& errorMessage, STring s_SRS, String t_SRS)

at Create.Manager.GDALManager.Convert(ImageFile imageFile, String& errorMessage)

GDALManager Code referenced:

try
{
    if (imageFile.CoordSys == CoordSys.UTM)
    {
        warpResult = GDAL.Services.GDALService.Warp(source, tempImageFileName, out fnErrorMessage, null, "NAD83");
    }
    else
    {
        warpResult = GDAL.Services.GDALService.Warp(source, tempImageFileName, out fnErrorMessage, imageFile.CoordSys.Name, "NAD83");
    }
}
catch (Exception ex)
{
    fnErrorMessage = "An error occurred while warping the file: " + ex.ToString();
}

Source and Destination Path: "\\DataServer\Images\Raster\Raster_Seamless"

All source and destination file paths are the same for both Dev and Production.


texts: Add value to entire column in attributes

texts:
Is there a way to add a value to the entire column? I don't think i can use raster calculator because it's a text i want to add, the letter "B".

Here is a picture.

This is what happends when i use field calculator.

I forgot the """""""". That's why it didn't work.


texts: Problem in python when creating a QGIS plugin

texts:
I have a problem with the command make in the python.

I used the:

Run->cmd and go to the location of your plugin builder

(i.e: c:/users/rafik_pc/.qgis2/python/plugins/pluginbuilder) type the following command:

pyrcc4 -py2 resources.qrc -o resources.py

(I used -py2 because I'm using python 2.7)

restart QGIS and DONE.

But the problem persists.

How can I resolve this?


texts: Graphic Layer graphics redraw incorrectly after map zoom event

texts:
I have a web map application which adds arrows onto street segments signifying direction as a Graphics Layer. The direction of each arrow is calculated by this process:
1) build an array of the street segments from a query
2) find the from and to node of each segment, and get the geometry of each
3) calculate the slope of the line from these two nodes
4) use the slope to calculate the angle of the arrow graphic, which is drawn on top of the line segment st its midpoint.

This works as expected after I added a short timeout to the function.  Before I added a timeout, it appeared that some arrows were drawn before each new angle was calculated, resulting in  some arrows being drawn in the wrong direction (although all arrow were in the right location).  Since I added the timeout, I do not have this problem on the initial graphics add.

The problem I have now is that if I zoom in or out of the map, many of the arrow are redrawn incorrectly.  I believe the same thing is happening here as was before I added the timeout; the arrows are drawn before the angle for each arrow is calculated.

How can I fix this? Can I somehow disable the redraw function, or do I need to rerun the query and angle calculation?

I am using ArcGIS Javascript API 3.8.


texts: How to identify to which polygon a point belongs to?

texts:
I have 3-4 polygons containing lots of lon/lat points. I want to identify to which polygon a specific point belongs.
Any help in this matter is appreciated.


texts: Different Packaging Error Problem for General function failure

texts:
I am using ArcMap 10.3 and while using Share As > Map Package, I was trying to save the package to disc. There were no errors when I used Analyze. So I hit the Share button and after a little bit of work, the “Save Map Document” pop up window came up. Immediately after saving it I got a “packaging error_general function failure”

First, here at StackExchange, a similar issue was resolved using the Check Geometry function. However, after checking both my layers' geometry, I didn't find any errors.

Second, searching online I found that another user fixed this problem when two feature classes were using the same data source. I deleted one of the layers that I wasn’t using to make sure the data was being drawing from only one table, but I'm still having the problem.

It seems this broad error can indicate a variety of problems.

How do I troubleshoot to figure out why I’m getting this error?

It would be easier if it was a geoprocessing tool creating the issue.


texts: Use point shapefile to spatial query multi polygon shapefiles with Python?

texts:
I have an address point layer and would like to spatially query this layer using various polygon shapefiles to intersect this information so that each address contains the attributes from these polygon layers.

Is this possible using non ESRI Python? I have used ArcPy but don't know if this kind of thing is possible with open source Python options? If so is the Python script fairly simple? Possible examples?


texts: Converting Google Maps Lat Long to Square meter area

texts:
I want to convert a set of 5 Lat Long (Directly taken from Google Maps) into Sq Meter area.

As 5 points are there I have found the formula on Converting Longitude and Latitude Coordinates to Square Miles?

area = rad(x2 - x1) * (2 + sin(rad(y1)) + sin(rad(y2))) + rad(x3 - x2) * (2 + sin(rad(y2)) +
+ sin(rad(y3))) + rad(x4 - x3) * (2 + sin(rad(y3)) + sin(rad(y4))) +
+ rad(x5 - x4) * (2 + sin(rad(y4)) + sin(rad(y5)))

area = abs(area * 6378137.0 * 6378137.0 / 2.0)

but this formula is not giving me the correct area.

I am putting the latitude in x & Longitude in y
also as i vary the digits after decimal places it changes the answer (obviously)
(ex. 20.011619, 73.762870 have 6 digits after decimal)

So what is the correct method to use this formula with lat, long? Do I want to convert lat-long into another measure to get the area correctly?


texts: Select and Delete from Multiple Layers at once in QGIS

texts:
I have imported a CAD drawing into QGIS as a GPK. It has imported everything that was in CAD. Now I want to delete features that I don't need. See below:

I know how to turn each layer on and then select the features in that layer and then delete them, but as there are many layers it will take a lot of time to delete it one by one. Is there a way to select all of this together and delete it?


texts: Make a new shapefile from parts of two shapefiles

texts:
I have tried all day to connect two shapefiles, crop a part of them, and save it into a new Shapefile.

In the image you can see two Shapefiles and inside of the square is the part which I want to export into a new shapefile.

It is possible?


texts: How GEE align pixels and projection when reduce resolution?

texts:
Here is the thing: I want to count how many Landsat pixels 30x30m were included in a 500x500m MODIS pixel. It is basically a zonal statistics, and I use reduceResolution to compute. But the issue is: in theory, there should be 500x500/30x30 = 278 landsat pixels in a MODIS pixel. But when I compute use count reducer in GEE, the number is 403 ~ 424, larger than the theoretical value. When implemented in arcgis, the sum zonal statistics returns 272 ~ 289, which makes more sense considering the actual MODIS pixel size is < 500 and Landsat pixel size is < 30 as well. So my question is how GEE handle the pixel resolution and alignment when implementing zonal statistics? Also, how can I get the value close to the theoretical 278? Is there a projection problem hidden in this issue? As GEE projection tutorial stated we should rarely concerned about projection when computing. I also tried reproject all layers to the same coordinates, but zonal statistics results are around 415.
here is my snippet using GEE example.
// Load a MODIS EVI image.
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A1').first())
    .select('EVI');

// Display the EVI image near La Honda, California.
Map.setCenter(-122.3616, 37.5331, 12);
Map.addLayer(modis, {min: 2000, max: 5000}, 'MODIS EVI');

// Get information about the MODIS projection.
var modisProjection = modis.projection(); // what's the difference betwen .projection() and crs()
print('MODIS projection:', modisProjection);

//var modisCrs = modis.crs(); // modis.crs is not a function
//print('MODIS crs', modisCrs);

// Load and display forest cover data at 30 meters resolution.
var forest = ee.Image('UMD/hansen/global_forest_change_2015')
    .select('treecover2000');

var forestCount = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.count(), // count gives value 403~424
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

// Display the aggregated, reprojected forest cover data.
Map.addLayer(forestCount, {max :80}, 'forest cover count at MODIS scale');

Here is screen shot of images when use ee.Reducer.count(). the black pixel value is around 403, and the white pixel value is around 424. Way large than the theoretical 278.


texts: Labeling only one duplicate attribute in QGIS without editing original data

texts:
I have several points (e.g. 4 bus stop positions) in one area. They all belong to the same station, so they all have the same name attribute. To avoid a mess of lables showing all the same name I only want to label one feature.

I dont want to create new (virtual)layers or columns in my file to keep it simple and universal to use; Just use them as they are. So an expression should do it.
Data is taken from osm (geofabrik).


texts: Programmatically modify MapTip Expression

texts:
If I load a layer via a shapefile and set

((ILayer)myLayer).ShowTips = true;

How can I change the tip content (I want to point it to a different column or have a expression of multiple columns) programatically with ArcEngine 10 and C#?


texts:
I like to implement a special WMS-layer with leaflet() in R. But I have no success with the WMS-layer I like to use.

library("leaflet")
library("sp")

# this example from the tutorial works
leaflet() %>% addTiles() %>% setView(-93.65, 42.028, zoom = 4) %>%
  addWMSTiles(
    "http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi",
    layers = "nexrad-n0r",
    options = WMSTileOptions(format = "image/png", transparent = TRUE)
)

# but for my implementation it don't work
leaflet() %>% addTiles() %>% setView(lat = 53.8, lng = 12.6, zoom = 8) %>%
  addWMSTiles("http://www.umweltkarten.mv-regierung.de/script/mv_a3_wasserschutz_wms.php?",
              layers = "t2_bib_p",
              options = WMSTileOptions(format = "image/png")
)

I assume that there is a problem with the Coordinate Refernce System. Because the WMS-layer I like to use, have not a standard CRS. So I try to expermented with this WMSTileOptions():

crs=CRS("+init=epsg:5650")

The idea is to use WMS-Layers as background images for my vector data. The vector data come mainly from shape-files and GPS-data (text files). The purpose is to create maps which i can use in a report.
If you know other possibles beside leaflet to use WMS in R there are also welcome?


texts: Sorting table in polygon shapefile in order of intersecting line?

texts:
I have 119 parcels that intersect a polyline.  I want to be able to sort the table to see if any of the land owners in the parcels have have consecutive ownership along the line.

Is there a way to do this, or some other process to see the order?


texts: Using Construct tools to construct a feature from line to polygon?

texts:
Is it possible for me to build a construct a line to polygon? For example, when you need a break or divide one polygon into two polygons. You had the line from a polygon to Line and you use that line to construct the feature that polygon to break two.

Currently, I tried to use this method using the Advanced Editing and when you selected two features: one line and one polygon even though the icon on the Advanced Editing shows no grayscale.

Is there a way around to do this ? It would save A Lot of time to use the "Cut Polygons" or any other tools.

Revised to explain more in detail.


texts: ArcMap function in field calculator not working

texts:
Trying to return everything except the last element in a field after splitting it. Using ArcMap 10.2.2 in the field calculator.

!ST_NAME!.split(" ")[:-1]

This returns a null value in attribute table. If a I use the same approach in python interactive window it seems to work:

a = "Johnny be good"

>>> b = a.split(" ")[0:-1]
>>> b
['Johnny', 'be']


texts: Can't assign a numeric value to a field that is set to

texts:
Within a feature class there's a number of records with a specific field (type double) set to . I've set up a cursor to assign those value a numeric number. However, I am getting an error.

This is the code:

import arcpy

working_folder = os.path.normpath(os.path.join(os.path.dirname(__file__), '..'))

#Input variables.
fc_Opt_All = os.path.join(working_folder, 'DATA\Optimity_Equipment.gdb', 'Optimity_All')

#Set initial number
last_location_id = 10000

field_names = ['LOCATION_ID']
i = 1
with arcpy.da.SearchCursor(fc_Opt_All, field_names) as cursor_1:
    for row_1 in cursor_1:
        if row_1[0] is None:
            row_1[0] = last_location_id + i
            i = i + 1

The error I get is the following:

TypeError: 'tuple' object does not support item assignment

The field 'LOCATION_ID' is numeric (double). I don't understand why it's saying it's a tuple. Any ideas?


texts: Getting pixel values from lat/long values using ENVI

texts:
I am using ENVI software to process satellite image. I have two images from two different satellites and I want to get the value of the same point in the two images. I got random points from the first one and I know the lat/long of each point.
How can I get the value of these points in the other image?


texts:
Im trying to produce a contour map using QGIS and exporting using Mapbox. The map is rotated 90 degrees in QGIS and it returns to its normal orientation if I use Leaflet. Unfortunately, using the Mapbox output, the transparent parts of the layer come out black.

Any suggestions as to how to fix this?


texts: Setting project QGIS variables with PyQGIS

texts:
In a plugin development I have to store in a project variable the plugin version.
To test if the variable already exists I wrote:

project = QgsProject.instance()
scope = QgsExpressionContextUtils.projectScope(project)
if (scope.hasVariable('plugin_version')):

To set the variable I use:

QgsExpressionContextUtils.setProjectVariable(project, "plugin_version", "0.1")

But if I try to add or set the variable through the scope (scope that I have to use to test if the variable exists) nothing happens.

scope.addVariable(QgsExpressionContextScope.StaticVariable("plugin_version", "0.1", True, True, "GisFIRE version"))
scope.setVariable("gis_fire_version", "0.11")

Why is this behaviour happening? How can I make a variable read-only?
QGIS version is 3.10


texts: How to merge pixel values to the nearest maximum value based on area in a classified raster?

texts:
We have a raster image with 9 discrete values. Generated the image using SAGA's Terrain Analysis -> Morphometry -> TPI classification technique. Now, like to merge pixels which occupy certain area (say less than 1 sq km) to the nearest maximum value. Eg, in the image - like to merge the white encircled ones to the nearest maximum value - 5 (Violet).

Because, if i try to vectorize this image, it takes too much time (and not getting finished) to convert those pixels into polygons as there are so many polygons with less than 1 sq km. Our aim is to get a vector file from the classified raster.

Is there any workaround like, split the raster into four and then do the same raster -> vector conversion ? But still, it will have too many polygons to deal with i guess


texts: I'm running R 3.4.4 and QGIS 2.14.9 on a linux server. Is there a possibility to utilize QGIS functions inside R - like the ArcGIS-R-bridge (which sadly only works on Windows)?


texts: How to get current date "YYYY-MM-DD" in string format in GEE script

texts:
I tried to get current "YYYY-MM-DD" in Google Earth Engine script.
I succeeded to get current "YYYY-MM-DD" in "Date" format in GEE script.
But I want to convert it to String format for setting as a parameter for "value:" of ui.Textbox() function.
Could you please help me ?
Code of verification of mine is follows.

var now = Date.now();       //getting time at UTC
var eeNow = ee.Date(now);   //converting number to object format

print("now:",typeof(now),now);
print("eeNow:",typeof(eeNow),eeNow);

var y=eeNow.get('year');
var m=eeNow.get('month');
var d=eeNow.get('day');

print("year:",typeof(y),y);
print("month:",typeof(m),m);
print("day:",typeof(d),d);

var Today01=ee.Date.fromYMD(
                eeNow.get('year'),
                eeNow.get('month'),
                eeNow.get('day'));
//                eeNow.get('day')).toString();

print("Today01:",typeof(Today01),Today01);

var Today02 = y.toString()+'-'+m.toString()+'-'+d.toString();
print("Today02:",typeof(Today02),Today02);

//
//  Final aim : to set "value:"parameter with current "YYYY-MM-DD" in String format
//
//var TxBoxt = ui.Textbox({
////////////////////////value: '2020-03-16',
//  value :Today01,
//  style: {width : '90px'},
//  onChange: function(text) {
//    var period_former_end = text;
//  }
//});


texts: Add a where condition (attribute) on a ArcGIS REST API FeatureService Layer applyEdits request

texts:
I am looking for a way to add a where condition (attribute) on a ArcGIS REST API FeatureService Layer applyEdits request

The idea is to filter the applyEdits executed on client side (ArcGIS JS Web App). The user can see everything and modify as well but the request must be intercepted and checked if the applyEdits are allowed against an attribute filter (where condition)
my environnment is based on ArcGIS Server 10.6 and ArcGIS JS 3.22.


texts: How to create a grid for a specific area and scale?

texts:
I have a municipal boundary with parcels that I need broken down into smaller maps to better view the parcels.  My client wants each smaller map at a scale of 1:10,000.  I plan to used data driven pages for this series of maps.  So I know I need to create a grid.  I want to use the Grid Index Features tool, my problem is that I'm unsure of what to use for my parameters. I don't have a specific number of columns or rows and I don't know how big of a polygon I want.

So how should I set up my parameters?  And is there a better way to go about doing this?

I have ArcMap 10.2.2 - Advanced


texts: Collapse dual carriageway lines in PostGIS

texts:
I have OpenStreetMap data, which represents larger roads (e.g. dual carriageways, aka divided highways) as parallel lines.

The question is answered for ArcGIS or Skeletron here, for example using ESRI's Merge Divided Roads or Collapse Dual Lines to Centerline tools.

Does anyone have suggestions for how to do this within PostGIS? I am working on a server without ArcGIS, and would like to avoid the overhead of transferring the data back and forth across machines.

Update
The linked question isn't really answered, and is looking for a solution that is not PostGIS. This should be kept open.


texts: Addition and subtraction of rasterbrick objects in R

texts:
I have two rasterbrick objects showing precipitation. Both have longitude and latitude dimensions of 144 x 72. The first (brick1) has 478 monthly time layers (39 years, 10 months) running from Jan 1979 to Oct 2018. The second (brick2) has 12 monthly time layers (Jan through Dec) representing the monthly mean values for 1981-2010. I want to calculate an anomaly by subtracting the geographically appropriate monthly mean from each data point. I did this:

anom <- brick1-brick2

R did not complain - but I am not sure if I got what I wanted.

My hope is that R will subtract the January through December means first against 1979, then it will recycle those means identically through each year until 2018, and in 2018 it will apply only the first 10 means (Jan through October). I have not been able to find any documentation that indicates whether this is the way it is actually done. I have seen some statements to the effect that some operations recycle smaller objects against the larger ones, but nothing specifically about rasterbricks. Does anybody know?


texts: Downloading from WFS - no data in returned XML

texts:
I want to download some information about lots using WFS hosted by http://geoportal.gov.pl/ .

I've tried download data with this link:

http://mapy.geoportal.gov.pl/wss/service/pub/guest/G2_GO_WFS/MapServer/WFSServer?service=WFS&version=1.1.0&REQUEST=GetFeature&typename=Dzialki&srsname=EPSG:2180&BBOX=30,30,30.001,30.001&

Unfortunately returned XML doesn't contain any interesting data...

<wfs:FeatureCollection xsi:schemaLocation='http://mapy.geoportal.gov.pl:80/wss/service/pub/guest/G2_GO_WFS/MapServer/WFSServer http://mapy.geoportal.gov.pl:80/wss/service/pub/guest/G2_GO_WFS/MapServer/WFSServer?request=DescribeFeatureType%26version=1.1.0%26typename=Dzialki http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd' xmlns:G2_GO_WFS='http://mapy.geoportal.gov.pl:80/wss/service/pub/guest/G2_GO_WFS/MapServer/WFSServer' xmlns:gml='http://www.opengis.net/gml' xmlns:wfs='http://www.opengis.net/wfs' xmlns:xlink='http://www.w3.org/1999/xlink' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance'>
<gml:boundedBy>
<gml:Envelope srsName='EPSG:2180'>
<gml:lowerCorner>133223.25099999999 171678.07000000001</gml:lowerCorner>
<gml:upperCorner>774923.77999999898 861902.51000000001</gml:upperCorner>
</gml:Envelope>
</gml:boundedBy>
</wfs:FeatureCollection>

I've tried variety of coordinates and BBOX values, but I've always got the same XML in response.

How can I get data from geoportal?


texts: Calculating the number of overlapping polygons in QGIS

texts:
I have several polygons drawn in Google Earth, saved as .kml files. They represent the same feature (a lava flow) as identified by different people on satellite images (so one polygon = one person). My goal is to build a map as this one, where the color represents the number of overlapping polygons, here from 0 to 6:

I've done this with Matlab, by "rasterizing" (using the inpolygon function and a meshgrid). The problem is that I lose the vector information. For instance, I can't calculate the area of the polygon where the 6 flows overlap (I can sum all the pixels with the same value, but this is imprecise because of the pixel size). Plus I get some distortion due to the shift in coordinate system.

I figured I'd better do it with QGIS but this is where I get stuck. I managed to load the polygons, and that's pretty much it... I think I should use a geoprocessing tool like intersect but it takes only two layers as input. Is there any way to do this?


texts: How do I offset a point created at the midpoint of a line in QGIS?

texts:
I am trying to accomplish something similar to this question, but in QGIS.
I am trying to create line segments from a polygon, while preserving on which side of the line, the polygon 'interior' is.
My idea is to create a point at the midpoint of each line segment, and then offset it onto each side of the line. I can then use a tool to check if the point is contained within the polygon, and mark each segment as "Interior_Left" and "Interior_Right"
Steps:

Polygon > line
Explode line to segments
Interpolate point at midpoint.  $length/2
Where I need assistance - Offset the point perpendicularly to both sides of the line, check if point is within polygon, and mark accordingly.

Now I am aware of GRASS plugin v.transform, but I don't know how I would offset it in relation to the azimuth of the line (perpendicular to the line segment)


texts: Finding the side of a square that a point intersects with R (sp)

texts:
I have a point and a square polygon. Both are sp classes.

The point intersects one side of the polygon, and I need to know which side of the square it intersects (i.e. north, east, south, or west). The square almost lines up with the compass, so there are clear n, e, s, w edges, but it doesn't line up perfectly.

The diagram illustrates this; in this case the point intersects the northern edge.

Any ideas how to work out which edge it intersects? On the surface it seems like a simple problem, but a solution has eluded me.


texts: Overwrite feature service error: 'object has no attribute'

texts:
I've published a feature service from a zipped shapefile. Now when I want to rewrite this service, I get this message.
AttributeError: 'str' object has no attribute 'type'

I don't know when there could be a problem? My item id is correct.


texts: Create Raster Surface from Points and Equations in ArcGIS

texts:
I want to use point data to create a raster surface which distributes the values over space, as defined by an equation, and NOT by interpolating.

I have an equation to describe how the characteristic changes  to describe a characteristic that decreases with distance from source e.g. y = x - 20 log (d) -11

Where d is the distance from the source, x is the value of the point data and yis the value for a point in the raster.

Thus the value of the raster decreases radially from each point. When values from two different points meet I would also like define how the values are added together. As they are not a simple addition of individual values.

I've included this graph to illustrate my question. Two point values (shown in red) a distance apart. The values the raster cells should contains are shown in blue. The shape of the blue line is very different to what would have been produced from an interpolation.

Is there a tool or set of tools that could be used in ArcGIS 10.2 to solve this problem?


texts: Shapefile not showing in proper location

texts:
Working in QGIS 3.0.1; I just switched to QGIS and am trying to get used to it.

I have a shapefile from SW Maps app of a track walked on hiking trails, and am attempting to make a basic map using this and a topo quadrangle for the area. However, while both are in same CRS (WGS84), the shapefile shows up in a completely different area. What am I doing wrong?


texts:
Can anyone explain, as simply as possible, how to calculate the visible sky percentage?

I have done a satellite survey and have worked out the height of buildings and the distance between them.

How do I use that data to calculate the visible sky percentage according to where I was standing?


texts: Which true type editor do you use to create/modify point symbols for GIS?

texts:
Which true type editor do you use to create/modify point symbols for GIS and what are the advantages and disadvantages of the tool?

I am looking for a true type editor to create/modify custom point symbols used in GIS. Here are some examples of symbols I want to create:

I am not familiar with the structure of a true type. But I guess that a true type contains many features which are not interesting for GIS (e.g. kerning). So the true type editor should be easy to use by a GIS analyst. I tried to get informations about true type editors by reading some descriptions in WWW (e.g. FontCreator). Before I test (or buy) a program I first want to get informations which editor is good for “GIS symbols”.

Both commercial and free solutions are welcome.


texts: Getting latitude and longitude raster from DEM using ArcGIS Desktop?

texts:
Could anyone provide with a way to calculate latitude and longitude from DEM, for each raster cell?

In previous posts about this, I have seen that some suggest calculating lat and long in certain sample point locations. This doesn't work for me, as I need complete rasters.

Also, I use ArcGIS 10.2 version so $$XMAP, $$YMAP does not work for me.

Is there any other way inside ArcGIS to get lat and long?


texts: compute average surface temperature from nc4 file with each lat long

texts:
I have nc4 file and I have one shapefile. So I want to read temperature for each latitude and longitude. I used that code-

import numpy as np
from netCDF4 import Dataset
import geopandas as gpd
import shapefile
data = Dataset("file.nc4","r")
//result -dimensions(sizes): lon(1440), lat(600), time(1), bnds(2)
          //variables(dimensions): float32 lat(lat), float32 lon(lon),float64 time(time), float64 time_bnds(time,bnds), float32 AvgSurfT_inst(time,lat,lon)

sf = shapefile.Reader('DISTRCT_BDY.shp')
for row in sf.shapeRecords():
    #print (row.shape.points)
    print (row.record)
    //result - ['AHMADNAGAR', 'MAHARASHTRA', 1355040.74938, 16554473772.2]

for temp in data.variables['AvgSurfT_inst'][:] :
        for tmp in temp:
            print (tmp)

so I am confused how to use lat and long in nc4 file to compute temperature.