texts: Performing Supervised Classification on Sentinel Image using ArcGIS Desktop?

texts:
Is it possible to perform supervised classification on Sentinel Satellite Image with the minimum error. I compared two raster, one being the manual classification of an area using Photoshop CS6 and the other being Supervised Classification with Interactive Supervised Classification (ArcGIS). The pixels that didn't match represented more than 70% (Confusion Matrix). So, this is the farthest I've gone. I need to know if there is more that I can do with only ArcGIS Desktop.

This is the Manual Classified Image with Photoshop CS6.

This is the raster using Interactive Supervised Classification (This tool uses training polygons)

As you see there is a big space that should be colored with violet in the second image. Obviously in the future using CS6 should not be what we always use but this is only used by compare and see if Interactive supervised classification works.


texts:
I'm working with the Leaflet JavaScript library. There should be an option to click on my map and view a popup. I already have a code for that:

var popup = L.popup();

function onMapClick(e) {
popup
    .setLatLng(e.latlng)
    .setContent(e.latlng.toString() + "
Visit Google")
    .openOn(map);
}

map.on('click', onMapClick);

This code is working, but my question is: how can i put an link to that popup? For example there should be  a link called "Visit Google". When you click on that link it should open www.google.com (surprise). I think it should be in the .setContent(e.latlng.toString()) tag, because I can add there a text like Visit Google (as you can see).

Any ideas ?


texts: Calculating median in grouped data using ArcMap

texts:
I got stuck when calculating the median age in grouped data (age ranges) for each feature (each feature is a single administrative unit). I have about 18,000 features.

Field names indicate age ranges: 0-14, 15-64, 65-100. In the "median" field I would like to calculate the median age for each feature.
I found a website https://atozmath.com/StatsG.aspx?q=1 where you can calculate such a median for a single feature, for example the first feature from the attribute table above.

The page easily calculated the median for the first feature in the image, and the page also wrote how it calculated it, but I don't know how to write a universal script for a "median" field to calculate the median for each feature.


texts:
I have 3 .nc files that were given to me without specifying which projection was used for the (lat, lon) coordinate pairs. Is there an attribute or a way of figuring it out given the latitudes and longitudes that are written in the file?

It's especially important to me to know since I need to obtain the closest grid points to a list of locations that are in epsg:4326 and compare the data stored in each.

I am using Python3 and GeoPandas to store the geographical information and plot.

I am new working with .ncand geographical data but it is my understanding different crs might change the (lat, lon) pair, so they cant be compared if they are not projected onto the same coordinate system.


texts: Is there a way of importing categorised layer styles into style-manager? I have created a set of styles for certain layers based on attribute data, but am having to go via properties>styles>load style> then save as default. It would be extremely useful to be able to load default style types based on layer naming like in the smart styles, but i believe this only works for styles with single symbols.


texts: Near analysis with mixed feature types in ArcObjects

texts:
I am trying to emulate the Near (analysis) ArcGIS tool in ArcObjects. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00080000001q000000

I want to use Mixed Feature Types. I have a point, polygon and line feature class and I want to use Near to compare all three of them with each other. This can be done one featureclass at a time like so (pseudocode logic):

Create IFeatureIndex
Store my feature class in IFeatureIndex.FeatureClass
Use IIndexQuery.NearestFeature for every IGeometry in a feature class to get distance to the nearest feature

The problem I am running into is that I can only store one feature class in IFeatureIndex. I need to store three different feature classes (points, lines and polygons) to be able to properly execute my comparison. How can I do this? ESRI's feature that I linked above has managed to do this but there is unfortunately no way to peek at the code behind the geoprocessing tool that I know of.

Edit: I have solved this. See the answer below.


texts: How to increase performance of tile caches if vendor parameter filters are needed?

texts:
My application uses vendor parameters to filter the data. It takes some time to generate all tiles on the server. Therefor I'm trying to increase the performance of the application.

I'm using geoserver as my tileserver. I looked up if caching would be a solution, but it does not support vendor parameters. The only solution I could think of when trying to cache tiles is to filter it on the client, by having a separate layer for each possibility. However this does not seem very doable.

I welcome any suggestions on how to improve performace.


texts: How are Modifications and Deletions represented when looking at enterprise geodatabase delta (add, delete) tables?

texts:
I'm trying to understand how to make sure a delete is a proper delete and not a modification, and vice versa. I'm using ArcGIS 10.2 and Enterprise geodatabase on MS SQL Server.

See below SQL:

select delete_table.*, add_table.* from
  db.schema.a175 add_table
  full join
  db.schema.d175 delete_table
  on
      delete_table.SDE_DELETES_ROW_ID = add_table.OBJECTID
      and add_table.sde_state_id = delete_table.sde_state_id
 where objectReference = 'Asset1'

and this returns something like

In this query, Asset1 is meant to have 2 objects related to it. But the user made two objects and realised one was wrong, deleted it and added another object. So now the Add table has three objects and the Delete table shows one.

But reading this output, I would have guessed objectid 8908810 is modified because there is also an entry in the Add table, and if this was a Deleted object, then the joined Add record would be all nulls.

Using the Add and Delete tables, how do I distinguish between modified and deleted objects?


texts: intersect a layer of polygons with itself, keeping the max attribute value in the result in ArcMap

texts:
I try to explain it as simple as possible.
I have a few overlapping polygons in my input layer, they have a rank attribute. Now i want to intersect them with each other creating a new layer polygons. These polygons should have the rank attribute of the highest ranking polygon of the input layer.

How can i do that in ArcMap.


texts:
I have a CSV layer of points which I have succesfully imported into QGIS. My assets are a bit like wireless access points and I cannot have any of the points on the same channel within 2km. I therefore want to be able to click on a point and it automatically draw a 2km circle around it, and then have the ability to click on 4 points at the same time to ensure the circles do not overlap, if it does then click off a point which will make the circle dissapear and I try another point.
So far I have tried adding a 'proposed channel' column in my CSV and reimporting the data every time I change a single point channel but this is taking forever... Does anybody have any better ideas?


texts: Should I Krige my contour data or use a different interpolation method?

texts:
I am using QGIS to calculate the terrain slope, starting from vector contour line data. Using linear interpolation like r.surf.contour gives me a reasonable looking DEM, but is not really useful if I need realistic slope data. I have seen a lot about Kriging point data, but it is not something I have tried before and also I have a large number of points – more than you would normally use for Kriging. Should I Krige the point data or is there a simpler way of getting a smooth slope?

 1. The DEM linearly interpolated from contour data

2. Slope of the DEM — the linear interpolation is clearly visible

3. Nodes from contour vectors — too many for Kriging?


texts: Going to the zoomed view of a particular 3Dbuilding using an URL in cesium API

texts: I have a cesium application which shows 3d buildings for an area. What I want is if I give some URL it should go to the zoomed view of a particular building. Say for ex: if I give http://localhost:8383/Cesium_3D_Buildings/Apps/Sandcastle/gallery/GeoJSON_Hubli_3D.html/kimscentrallibrary it should go to the Kims central library building with the zoomed view. I know it can be done somehow by fly to but is is possible using a URL in cesium API?


texts: Unable to retrieve attributes after assignment

texts:
Long story short I cannot get the attribute(s) back for labeling after assignment. I must not be doing it right. I format the attributes as follows:

{'fipscode': '09003','state': 'ny','county': 'kings county','speccode': 'AWSB1U','callsign': 'WQGA906','fcc': 'Cellco Partnership','lease': 'NULL','leaseid': 'NULL','current': 'Cellco Partnership','future': 'Cellco Partnership'}

and assign them as follows:

feature.attributes = attrsList;

and when I attempt to retrieve the county field, or any field, it is undefined.

More Code Follows, Ideas?

function initMap() {

var wktReader = new ol.format.WKT();

var featureCollection = [];

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

    var feature = wktReader.readFeature(geometries[i]);
    var attrs = attributes[i];
    var textMsg = "";
    var state = "";
    var county = "";

    var mystyle = new ol.style.Style({
        fill: new ol.style.Fill({ color: 'rgba(0, 0, 255, 0.1)' }),
        stroke: new ol.style.Stroke({ color: 'blue', width: 3 }),
        text: new ol.style.Text({
            textAlign: 'center',
            textBaseline: 'middle',
            font: '12px Verdana',
            text: GetAttributes('fipscode'),
            fill: new ol.style.Fill({color: 'black'}),
            stroke: new ol.style.Stroke({color: 'white', width: 0.5})
        })
    });

    if (feature.getGeometry().getType() == 'Polygon') {

        feature.getGeometry().transform('EPSG:4326', 'EPSG:3857');

        // Assign the attributes
        var attrsList = GetAttributes(attrs);

        feature.attributes = attrsList;
        //feature.features = attrsList;

        //var looky = feature.get('fipscode');
        //feature.text = feature.get('fipscode');
        feature.setStyle(mystyle);

        //feature.setStyle(styleFunction(feature));
        // styleFunction()

        featureCollection.push(feature);
    }
    else if (feature.getGeometry().getType() == 'LineString') {
        feature.setStyle(new ol.style.Style({
            stroke: new ol.style.Stroke({
                color: 'red',
                width: 3
            })
        }));

        featureCollection.push(feature);
    }
    else if (feature.getGeometry().getType() == 'Point') {

        feature.setStyle(new ol.style.Style({
            image: new ol.style.Icon(/** @type {olx.style.IconOptions} */({
                anchor: [0.5, 46],
                anchorXUnits: 'fraction',
                anchorYUnits: 'pixels',
                opacity: 0.75,
                src: 'Icons/marker.png'
            }))
        }));

        featureCollection.push(feature);
    }
}


texts: Why is ArcPy script looping through bands slow?

texts:
I am working on a single tiff with over 20k bands, each representing daily temperature, and I'll be extracting and preforming simple raster math on roughly 3500 of the bands. The problem I'm running into is that each band is taking ~2 min to process which would equate to a 4 day running time. I am sure there must be a more efficient manner for processing this.

*edit: When running for 1 band in range (1, 2) I am still getting a ~2min process time.

*edit 2: Relevant information of file.
Columns, Rows:  179, 195
Number_of_Bands:  20454
Cell_Size__X.Y:  0.0625, 0.0625
Format:  TIFF
Pixel_Type:  floating point
Pixel_Depth:  32 Bit

I have run the same script on a subset with ~200 bands from the file above and the processing time per-band is nearly 2 seconds each.

>>import took 0:00:20.362000
>>set up environment took 0:00:08.220000
>>a full band loop took 0:01:43.342000
>>a full band loop took 0:01:42.302000
>>a full band loop took 0:01:41.980000
>>a full band loop took 0:01:41.359000
>>a full band loop took 0:01:41.072000
>>one season loop took 0:08:32.306000

import datetime
t1 = datetime.datetime.now()

import arcpy, os, arcinfo
from arcpy.sa import *

t2 = datetime.datetime.now()
print "import took %s" %  ( t2-t1)
arcpy.env.overwriteOutput = True

# Set the current workspace
arcpy.env.workspace = "C:\\raster\\folder"
arcpy.CheckOutExtension('Spatial')

# Get and print a list from the workspace
rasters = arcpy.ListRasters('*', 'TIF')
cutoffTemp = float(277.2)

t3 = datetime.datetime.now()
print "set up environment took %s" %  ( t3-t2)

for raster in rasters:
    t4 = datetime.datetime.now()
    print(raster)
    d = arcpy.Describe(raster)
    number_of_bands = d.bandCount #number of bands from Hist2005
    print(number_of_bands)

    arcpy.AddMessage("Processing raster {} for {} bands".format(raster, number_of_bands))
    listcomp90 = [] #create a list of output tifs, empty for each input raster

    for band_number in range(1, 6): #small range for test
        t6 = datetime.datetime.now()
        listcomp91.append(band_number)
        Hist = arcpy.Raster("{}\\Band_{}".format(raster, band_number)) #gets inidiviual bands for input
        arcpy.AddMessage('in the band number cycle and the current band number is {}'.format( band_number))
        output = arcpy.sa.Con(Hist >= cutoffTemp, sapHist) #if band values >= cuttoffTemp then temp
        Rasname, Rasext = os.path.splitext(raster)
        output.save(os.path.join(arcpy.env.workspace, "New{}_91_comp_{}.tif".format(Rasname, band_number)))
        listcomp91.append(os.path.join(arcpy.env.workspace, "New{}_91_comp_{}.tif".format(Rasname, band_number)))
        t7 = datetime.datetime.now()
        print "a full band loop took %s" % (t7-t6)

t5 = datetime.datetime.now()
print "one season loop took %s" % (t5-t4)


texts: Problem in displaying tile images using OpenLayers 3

texts:
I have problem with displaying my TMS layer using OpenLayers 3 (OL3).
This is my first time working with tile images this way so maybe it doesn't have to do anything with OL3, maybe I just don't understand it right.

However, I've downloaded sample .tif image from NaturalEarthData http://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/ and then I generated my map tiles (for zoom levels 0-5 ) using gdal2tiles command. Tile images are placed in nested folders following the zoomLevel/xTile/yTile.png schema.

I created simple map using OL3 and added tile layer to it. At first zoom level it's displayed fine.

However, when I try to zoom in it doesn't work good because tile images are not merged correctly.

Can you please tell me where I'm going wrong, I tried to follow this example http://tileserver.maptiler.com/#nasa/ol3 but I just can't achieve same thing with my data.

You can see my code below:

var extent = ol.proj.transform([-180,-85.051129,179.976804,85.051129],
  'EPSG:4326', 'EPSG:3857');
var center = ol.proj.transform([-0.011598000000006436, 0],
  'EPSG:4326', 'EPSG:3857');

var map = new ol.Map({
  layers: [
    // new ol.layer.Tile({
    //   source: new ol.source.MapQuest({layer: 'osm'})
    // }),
    new ol.layer.Tile({
      source: new ol.source.XYZ({
        url: 'static/website/NE1_HR_LC/NE1_HR_LC/{z}/{x}/{y}.png',
        extent: extent,
        minZoom: 0,
        maxZoom: 5,
        wrapx: false
      })
    })
  ],
  renderer: 'canvas',
  target: 'map',
  view: new ol.View({
    projection: 'EPSG:3857',
    center: center,
    zoom: 1
  })
});
map.getView().fitExtent(extent, map.getSize());


texts: Reading Elevation Data from a DSM file using Python

texts:
I am new to GIS processing and I have a DSM file with a .tif extension. I would like to read the elevation data from this file using python or QGIS.

I can view this file on the delair web application like below.

The objective is to identify the mounds you can see in the center of the image using elevation data , However, I am not even able to read it.
Things which I have tried:

Installed QGIS software and opened the DSM file in it and it appears to be a GUI  for the same thing as this website does , Not sure how to get this from an image to values.

Tried using gdal python package to open this DSM .tif file and it
throws exceptions.

Note: I am not from this field and have limited knowledge mostly only python, I am covering for a fellow employee who hurt his leg and wont be in for a couple of months and we really need to get started on this project.

If my total approach is wrong kindly let me know.


texts: Using multiple layers as point features in spatial join using ArcMap

texts:
I have 15 point feature layers and one polygon layer. I need to count the total number of points in the polygon using a spatial join, to then classify the different polygons based on the join count.

Is there any way to do this with another task, because the spatial join wizard only allows me to use one point layer?


texts: How to open and plot Globcolour chlorophyll image in R?

texts:
Level 3 globcolour chlorophyll images (download from hermes) are in netCDF format. Those files can be imported via ncdf package to R 3.1.

ch.nc = open.ncdf("E:/Data/globclour/L3b_20140101-20140131__GLOB_4_AV-MOD_CHL1_MO_00.nc")

How can I plot CHL1_mean? I follow guide by Mark but not getting any success.


texts: How to publish OpenStreetMap (OSM) layer on GeoServer directly?

texts:
Is it possible to directly publish OpenStreetMap (OSM) format data to GeoServer?

Can any one guide me to publish i.e. without loading the OSM data into MySQL or PostgreSQL.

Just publish the OSM format data from the folder where the file is located.


texts:
I see a number of posts about changing the resolution of an output of a map; however, my problem is that the map itself is lacking in adequate resolution. When working in the main QGIS program, I can zoom in and see much detail in my underlying aerial photo. However, when I make a map using Print Composer, I cannot zoom in without losing too much resolution. It becomes pixillated.

I have tried increasing the size of the map, but that doesn't seem to make any difference. My properties are set to "render", not "cache". Under "Variables\Composition\layout_dpi" it is 300; but I can't seem to change that, or don't know how.


texts: Simplifying stream network based on grid using QGIS?

texts:
I have a river network (blue line in the below image) and I need to simplify this based on a grid. I basically need to snap the river to the closest vertex or middle-of-side of the grid cell's bounds (grey dots) so the simplified river is that shown in red.

The ultimate goal is so that I can create this type of representation of the river in data (to input to a model), that describes the network as a series of river reaches (within grid cells) that are linked by pointing to their inflow grid cells.

"GridCell_1_1": {
    "RiverReach_1_1_1": { ... }
},
"GridCell_1_2": {
    "RiverReach_1_2_1": { "inflow": "RiverReach_1_1_1" }
},
"GridCell_2_1": {
    "RiverReach_2_1_1": { ... }
},
"GridCell_2_2": {
    "RiverReach_2_2_1": { "inflow": "RiverReach_1_2_1" },
    "RiverReach_2_2_2": { "inflow": "RiverReach_2_1_1" }
}

How can I achieve this using QGIS?


texts: arcpy MakeFeatureLayer - one workspace question, one output question

texts:
I'm trying to run a very simple python script to create a feature layer from a shapefile.

The original script read as follows (I added the two print statements at the end to verify upon completion that:

1) the output I was trying to create was created (Exists(the layer I was trying to create)) and

2) that the script was recognizing the path name by asking if a known folder existed

# Import system modules
import arcpy
from arcpy import env

# Set overwrite option
arcpy.env.overwriteOutput = True

# Set workspace
# arcpy.env.workspace = r"H:\Python\Python_Cookbook\data\TravisCounty"

# Make feature layer to hold selection.
Bldgs = r"H:\Python\Python_Cookbook\data\TravisCounty\BuildingPermits.shp"
Bldg_Lyr = r"H:\Python\Python_Cookbook\data\TravisCounty\Building.lyr"
arcpy.MakeFeatureLayer_management(Bldgs, Bldg_Lyr)

print arcpy.Exists(
    r"H:\Python\Python_Cookbook\data\TravisCounty\Building.lyr")

print arcpy.Exists(
    r"H:\Python\Python_Cookbook\data\TravisCounty")

The result from this first script was that no layer appeared in Catalog, and the print statements returned False and True. So the script recognized the existing folder, but did not find the layer it was supposed to create.

My second attempt was to comment out the workspace and to assign variables to the full paths instead. (see below).

# Import system modules
import arcpy
from arcpy import env

# Set overwrite option
arcpy.env.overwriteOutput = True

# Set workspace
# arcpy.env.workspace = r"H:\Python\Python_Cookbook\data\TravisCounty"

# Make feature layer to hold selection.
Bldgs =     r"H:\Python\Python_Cookbook\data\TravisCounty\BuildingPermits.shp"
Bldg_Lyr = r"H:\Python\Python_Cookbook\data\TravisCounty\Building.lyr"
arcpy.MakeFeatureLayer_management(Bldgs, Bldg_Lyr)

print arcpy.Exists(
r"H:\Python\Python_Cookbook\data\TravisCounty\Building.lyr")

print arcpy.Exists(
r"H:\Python\Python_Cookbook\data\TravisCounty")

This seemed to work better: My print statements returned True and True this time around.

HOWEVER: When I opened the ArcCatalog, I still couldn't find the new layer that was supposed to have been created.
Copying the full path (that just returned 'True' to the exists() call) into the Location bar of ArcCatalog returned a message that it was "...and invalid or non-existent directory."

ps. Is there a way to apply the codewrapper without manually entering four spaces before each line of code?


texts: Are there known issues about UPPER, LOWER SQL functions in ArcGIS 10, file GDB?

texts:
Platform: ArcGIS 10 SP3, .NET

I use:
IQueryFilter queryFilter = new QueryFilterClass();
with ArcGIS .NET SDK and build a query...

Does anybody had problems using UPPER(stringField) and LOWER(stringField) functions on esri file geodatabases? If so, please share links to some "confirmed-will-fix-in-2-years" pages...

I know we have some "old" warnings in our codebase about that. Like we though UPPER(stringField1)=UPPER(stringField2) doesn't work on file gdb. Last week I ran into it again but it seems I could not get even UPPER(stringField)=('STRING-IN-mixedCACE'.toUpper()) query to find anything in gdb.
UPPER(stringField)='STRING-IN-UPPERCACE' works, but the UPPER should be better named NOOP.

Example. There is a record in fgdb where stringField='Abc'. The query UPPER(stringField)='ABC' returns 0 records. (Error: Should be 1!)

PS. Such UPPER/LOWER queries seem to work correctly in ArcMap when using feature layer definition query.

Example 2.

Table (id, stringField)
1, 'Abc'
2, 'ABC'

QueryFilter: where UPPER(stringField)='ABC'
returns: (2, 'ABC')
should-return: (1, 'Abc'), (2, 'ABC')


texts:
I would like to do a Dijkstra routing using pgrouting and openlayers. I have the downloaded osm data. I am following the workshop given in this link. But since I am very new to this I found it difficult to understand this page.

Can anyone please explain me this page. Or provide me with some suggestions to do it in any other way.


texts: Display Lat/Lng coordinates from Oracle DB with ArcGIS HTML5 component

texts:
We have some experience in using ArcMap and ArcSDE, but this dates back to 2008, so we're not up to date with the recent developments in ESRI products.

Now we're dealing with a database that has some Lat/Lng coordinates stored as plain FLOAT columns. The client wants an HTML5 map that displays some information based on these Lat/Lng coordinates: some point annotations, as well as some arrows between pairs of Lat/Lng coordinates.

I proposed an architecture using Google Maps Javascript API v3 for the map component. I'd write a small web service that would deliver the Lat/Lng coordinates as JSON. In the browser, some small javascripts would display the Lat/Lng data in form of annotations and arrows. No big deal.

However, due to high-level decisions by the customer's management, the map component must be implemented using ESRI products, so my idea with Google Maps and the simple JSON service won't work as intended.

I see the following options:

Basic scripting: We're using roughly the same architecture as we would have used with Google Maps. So we need an HTML5 map component by ESRI that has a Javascript API, so we can display annotations and arrows via Javascript.
WFS service: We install/implement some kind of (OGC-compliant) WFS web service that reads the Lat/Lng values from the Oracle database and serves them as an OGC-compliant feature layer. Then, we're using an HTML5 map component by ESRI to consume that WFS service. No idea how I would implement the arrows between two Lat/Lng pairs.
ESRI service: Maybe ESRI provides a web service component with a generic data source, so we could wire that web service component to the Lat/Lng columns from the DB. This is basically the same architecture as in option 2, but with an ESRI service instead of some generic WFS component.

Thanks for any hints. We're new to the HTML5 world of ESRI, so we're just looking for some what's-possible and how-to-get-started hints.


texts: Fast(est) way to detect if point is in range of another point? (Geographic Coordinates)

texts:
I'm looking for a way to detect if 2 points are "near enough" each other for a match, and hopefully, if the range could be specified in metric system (each item may have different ranges, but there would be no change once it has been inserted).

I have already used Distance before for things like this, though the amount of data was significantly less, now I'm afraid it would be too slow for what I need.

Using simple logic, I thought about the fastest way to calculate that would be to do half the calculation already on the table, since one of the points has a constant range and location (neither can be changed later, worst case scenario, range can be changed, but then only need to update circle area), I could instead save it as a Circle format.

I'm not really sure, But as I looked more on it I reached a conclusion that usual GIS functions use the Cartesian Coordinate System for its calculations, which could work fine under a certain limited area, but because I would be working with GPS coordinates, I'm afraid that the farther away from points 0°,0°, the bigger the error margin may be too big (The radius could be from a few meters to a few Kms) and but I'm not sure if there are functions for Geographic Coordinates, or if anyone knows a way to draw a "deformed circle" accordingly based on the base latitude point.


texts: Converting BMP image into TAB MapInfo format using MapXtreme?

texts:
I am developing a plugin for MapInfo and I have BMP image with all the prop and I have rectbounds. I have done through GDAL using the below solution. But I am not getting proper output which I am expecting. I am getting some extra white patches in my TAB tab after converstion. I believe this is happening because of GDAL.

This is my code:

            Bitmap bmp = CreateChartTIF(path);

            GdalConfiguration.TryConfiguration();

            String wktProj = null;
            String tmpPath = @"E:\@@MapInfo_1\tmp.bmp";
            int a = Convert.ToInt16(rectBounds.Left);
            int b = Convert.ToInt16(rectBounds.Top);

            bmp.Save(tmpPath, ImageFormat.Bmp);
            System.IO.MemoryStream mem = new System.IO.MemoryStream();
            mem.Flush();

            String[] options = null;
            Gdal.AllRegister();
            OSGeo.GDAL.Driver srcDrv = Gdal.GetDriverByName("GTiff");
            Dataset srcDs = Gdal.Open(tmpPath, Access.GA_ReadOnly);
            Dataset dstDs = srcDrv.CreateCopy(@"E:\@@MapInfo_1\Output\2-0.tif", srcDs, 0, options, null, null);

            Osr.GetWellKnownGeogCSAsWKT("WGS84", out wktProj);
            dstDs.SetProjection(wktProj);

            double minx = rectBounds.Left;
            double miny = rectBounds.Bottom;

            double maxx = rectBounds.Right;
            double maxy = rectBounds.Top;

            double mapWidth = Math.Abs(minx - maxx);
            double mapHeight = Math.Abs(maxy - miny);

            double pixel1 = mapWidth / bmp.Width;
            double pixel2 = mapHeight / bmp.Height;

            //Top Left x, w-r pixel resln, 0, Top Letf y, 0, n-s pixel resln
            double[] geoTransfo = new double[] { minx, pixel1, 0, maxy, 0, pixel2 * (-1) }

            dstDs.SetGeoTransform(geoTransfo);

            dstDs.FlushCache();
            dstDs.Dispose();
            srcDs.Dispose();
            srcDrv.Dispose();
            bmp.Dispose();

            RasterProcessing.Convert(@"E:\@@MapInfo_1\Output\2-0.tif", "E:\\@@MapInfo_1\\Jonathan_Test.mrr", "MI_MRR");

Please find the screenshot below :

Any idea or any sample code? Is there any way to remove GDAL dependency and do this conversion using MapXtreme?


texts: I started making a vector layer from raster data but I see slivers (gap area) between my parks and roads layer. How can I overcome this problem? I planned to make a map with scale 1:5000.


texts:
I am at a complete loss on what I am doing wrong here.

I am trying to build a variogram.  No matter what i do i get errors.

My Data:

     id       lon      lat                  datetime         m1      m2       m3
1      1 -121.3501 41.86009 2016-03-01 00:00:00-08:00    43.200  3993.071  3993.071
2      2 -121.3412 41.94472 2016-03-01 00:00:00-08:00    38.800  4036.954  4036.954
3      3 -121.5138 41.99801 2016-03-01 00:00:00-08:00    24.600  4023.301  4023.301

After converting to a Spatial DataFrame

class       : SpatialPointsDataFrame
features    : 3196
extent      : -124.2696, -116.1304, 33.11097, 42.01195  (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables   : 5
names       : id,                  datetime,      m1,      m2,        m3
min values  :   1, 2016-03-01 00:00:00-08:00,     0.000,     0.020,  0.0003920282
max values  :   3196, 2016-03-01 00:00:00-08:00, 14360.500, -7470.500, -3.4202715765

This is just a small piece of the data. About 3000 rows in total

My dataframe setup

df <- read.table("input.csv", header=TRUE, sep=",")
df = na.omit(df)
coordinates(df) = ~lon+lat

Some samples of what I have tried and the error messages:

v = variogram(m3~1, data=df)

Error in variogram(m3 ~ 1, data = df) : unused argument (data = df)

v = variogram(df$m3~1, data=df, dX=0)

Error in variogram(df$m3 ~ 1, data = df, dX = 0) :
  unused arguments (data = df, dX = 0)

Since it keeps telling me these parameters are unused...I tried this:

v = variogram(df$m3~1)

Error in variogram(df$m3 ~ 1) :
  argument "X" is missing, with no default

All I am trying to do at this point is create a basic variogram on the m2 or m3 column of data.  m3 column is the m2 column with a normal score transformation done on it.  I just left the code out for that.

What am I doing wrong?


texts: Output GeoTIFF file is wrongly rotated to 90 degrees/ -90 degrees using xarray and rioxarray with NotGeoreferencedWarning warning

texts:
I am using this file. I wrote this code.
import rioxarray
import xarray

xds = xarray.open_dataset("LIS_HIST_201401010000.d01 copy.nc")

floodFraction = xds.FloodedFrac_tavg.transpose('east_west','north_south')
floodFraction.rio.set_spatial_dims(x_dim="north_south", y_dim="east_west", inplace=True)
floodFraction.rio.write_crs("EPSG:4326", inplace=True)
floodFraction.rio.to_raster('floodFraction.tif')

After running this script, my terminal shows this
NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific.
  s = writer(path, mode, driver=driver,

The obtained tif file looks like this


texts: OpenLayers: cannot register custom projections

texts:
I'm trying to load an OpenLayers map with UTM projection. If I try to register custom projection, the map becomes blank.

proj4.defs('EPSG:32640', '+proj=utm +zone=40 +datum=WGS84 +units=m +no_defs');

let options = {};
options.code = 'EPSG:32640';
options.units = 'm';

options.extent = [
  -98570.85212537996,
  2468701.5790765425,
  683268.1076887846,
  2874585.9453238174,
];

options.worldExtent = [];

options.axisOrientation = 'enu';

options.global = false;

//options.metersPerUnit = projectionOptions.metersPerUnit;
options.getPointResolution = ()=>{};

let projection = new ol.proj.Projection(options);
//ol.proj.addProjection(projection); // if I uncomment this the map becomes blank
//let ext=ol.proj.get('EPSG:32640').getExtent();
debugger;
//let center= ol.extent.getCenter(ext);
var olview = new ol.View({
  center: [312234.8270497762, 2636116.2213047906],
  extent:options.extent,
  projection: 'EPSG:32640',
  zoom: 8,
});
var map = new ol.Map({
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM(),
    }),
  ],
  target: 'map',
  view: olview,
});

Plunker

ol.proj.addProjection(projection); // if I uncomment this line, map becomes blank


texts: Is there a way to retrieve the shape of a styled object from a database?

texts:
I am working with geoserver on an oracle database.

I am trying to generalize my SLDs and make them more dynamic. The idea is to store all styling information in the database and let the SLD retrieve it dynamically.
I found out that this works well for fill, stroke and some other atrributes like opacity, size and rotation. As can be seen below where SLD_FILL and SLD_SIZE are the columns where the actual fill and size styling data is stored.

<PointSymbolizer>
<Graphic>
     <Mark>
     <WellKnownName>circle</WellKnownName>
         <Fill>
         <CssParameter name="fill">
         <ogc:PropertyName>SLD_FILL</ogc:PropertyName>
         </CssParameter>
         </Fill>
     </Mark>
     <Size>
     <ogc:PropertyName>SLD_SIZE</ogc:PropertyName>
     </Size>
</Graphic>
</PointSymbolizer>

I have tried this also with the shape of a marker where the column SLD-SHAPE stores the text circle:

<PointSymbolizer>
  <Graphic>
     <Mark>
     <WellKnownName>
     <ogc:PropertyName>SLD_SHAPE</ogc:PropertyName>
     </WellKnownName>
         <Fill>
         <CssParameter name="fill">
         <ogc:PropertyName>SLD_FILL</ogc:PropertyName>
         </CssParameter>
         </Fill>
     </Mark>
     <Size>
     <ogc:PropertyName>SLD_SIZE</ogc:PropertyName>
     </Size>
  </Graphic>
</PointSymbolizer>

This does not work. Does anyone know if this is possible?


texts:
I have ten shapefiles containing polygon features.

I would like to produce a single polygon layer containing portions of the input polygons where no fewer than two of the aforementioned feature classes intersect.  Any two (or more) polygons which intersect will have those portions included in the output feature class.

Additionally, each of the input feature classes contains a unique field name representing itself.  I would like for the output feature class to contain EACH of the unique field names from the input feature classes, as well as the value contained in the attribute, or 0 if the polygon did not intersect.

So, if a polygon from feature class A and feature class B intersect, that intersecting portion will exist in the output feature class.  8 of the 10 attributes in the output will contain the value "0", and the attributes for the fields representing feature classes A and B will contain the values from those input features.

I hope that makes sense.

I would like to accomplish this task using ArcGIS.


texts:
Is there anyway in QGIS that I can select all the outermost features on a map? On an island that would be all that have a coastline.

edit: This is what I'm trying to select:


texts: How are labels dynamically placed on a line or polygon?

texts:
I was wondering if anyone knew the math involved, or resources for learning the math involved, in finding the ideal position, scale, and rotation of a label for a line or polygon on a map?

The most I've got so far is finding the centroid for a polygon (though it's not ideal for concave polygons).


texts: QGIS Random Points Inside Polygons Based on Attbribute Value?

texts:
How can I set Random Points Inside Polygons so that the number of points created inside each feature is based on an attribute value for that feature?

Based on other answers I see elsewhere on Stack Exchange, there should be an option for this in the GUI but I'm not seeing it.

I'm using QGIS version 3.10.6 (A Coruña).


texts: Calling GRASS functions from external Python?

texts:
I have GRASS 6.4.3 svn installed on my machine running Windows 7. Also in my machine, I installed Python 2.7.3 for Windows. Now I need to call GRASS libraries from Python and not within the Python shell in GRASS GUI itself. I set the environment variables in the env.bat for this. But connection fails. This is the env.bat I have:

set GISBASE=C:\Program Files\GRASS GIS 6.4.3svn

set HOME=%USERPROFILE%

set LD_LIBRARY_PATH = %GISBASE%\lib

set PATH= %GISBASE%\etc;%PATH%
set PATH=%GISBASE%\Python27;%PATH%
set PATH=%GISBASE%\lib;%GISBASE%\bin;%GISBASE%\extralib;%GISBASE%\msys\bin;%PATH%
set PATH=C:\Python27;%PATH%

set PYTHONLIB= C:\Python27;%PYTHONLIB%
set PYTHONPATH= %GISBASE%\etc\python;%PYTHONPATH%

rem Environmental variables for GRASS stand-alone installer

set GRASS_WISH=%GISBASE%\extrabin\wish.exe
set GRASS_PYTHON=python

set GRASS_PROJSHARE=%GISBASE%\proj
set GRASS_HTML_BROWSER=explorer
set GRASS_SH=%GISBASE%\msys\bin\sh.exe

set PYTHONHOME=%GISBASE%\Python27
set GDAL_DATA=%GISBASE%\share\gdal
set PROJ_LIB=%GISBASE%\proj
set GEOTIFF_CSV=%GISBASE%\share\epsg_csv

set PATH=%GISBASE%\msys\bin;%PATH%
set PATH=%GISBASE%\extrabin;%GISBASE%\extralib;%PATH%
set PATH=%GISBASE%\tcl-tk\bin;%GISBASE%\sqlite\bin;%GISBASE%\gpsbabel;%PATH%
set PATH=%GISBASE%\bin;%PATH%

rem Path to the python directory --
rem set PYTHONHOME=%GISBASE%\Python27
rem if "x%GRASS_PYTHON%" == "x" set GRASS_PYTHON=python

Can somebody help me in rectifying the problem?


texts: Polygon creation fails: "The feature cannot be added because it's geometry collapsed due to intersection avoidance"

texts:
From one moment to the other I cannot create any polygon anywhere on the map, independent from the layer I am in, even if I create a new shape-file. I have checked all snapping options, enabled them, disabled them completely, restarted everything, etc. but no matter what I do this message pops up:

  The feature cannot be added because it's geometry collapsed due to intersection avoidance

It's super weird because even if I create it in absolutely blank space it tells me that the geometry collapsed. It occured when I had print composer open and modified a polygon. I moved the knots of one polygon and then wanted to create a new one which did not work. I checked the polygon of which I had moved the knots 400 times...there is no flaw and no intersection. Frustrating.

I am using QGIS 3.4 Madeira


texts: GRASS plugin Error: Unable to properly access 'C:\...'

texts:
I received this messsage in the log when trying to run the GRASS plugin for various functions:

I am using QGIS 3.10.9 with Grass7. The message warning says:
"WARNING: Concurrent mapset locking is not supported on Windows
ERROR: Unable to properly access 'C:\Users\Alex\AppData\Local\Temp\processing_34d0f624b1b2411eb5d7d59c5b3f9aad\grassdata\temp_location\PERMANENT.gislock'."
How would I fix this issue?


texts: Registering a SDE Service in 10.1 with Command Line

texts:
I am installing ArcSDE on a new development box.  I have executed the setup.exe, created an Enterprise database called "GDB123" and created a new sde service called "sde_123".  I need some clarification on how to register the service using the administrative commands.  The install guide says if your are using dbo-schema register the service with this command

sdeservice -o register -d SQLSERVER,MYSERVER\SSINSTANCE -r SDE_DBA_USER -v sa -i sdesvc

I do not understand what value I am suppose to type in for the -r and -v options.  I think -r needs to be the database name "GDB123" but I have no idea what to put for the -v option.  Heres what I have so far

sdeservice -o register -d SQLSERVER -r ?? -v ?? -i sde_123

How do I make this work?

ArcSDE 10.1 install for SQL Server on Windows 2008 RS box

Install Guid: http://resources.arcgis.com/en/help/install-guides/arcsde-application-server-microsoft-sql-server/10.1/index.html#/Creating_an_ArcSDE_service_on_Windows_using_the_sdeservice_command/00950000003v000000/

SDE Service Commands Options: http://resources.arcgis.com/en/help/install-guides/arcsde-application-server-microsoft-sql-server/10.1/index.html#/The_sdeservice_command/00950000004m000000/


texts: Can somebody help me in getting a random point which lies inside multipolygon,i.e I have a multipolygon, then i want to get a random point which lies inside that multipolygon


texts:
I'm trying to buffer points and scale the buffer around the points themselves. According to the documentation this can be done by supplying the factor as a point to the second argument and the origin to the third.

But my scaled buffers are not coming out correctly and have no extent according to QGIS. How can I scale the buffers around the points' origins?

Here's an example of what I have currently:

DROP TABLE IF EXISTS zones;
CREATE TABLE zones AS (
WITH
coords AS (
SELECT *
FROM ( VALUES
    (1, 101, 201),
    (2, 104, 204),
    (3, 109, 209)
   ) p(id, x, y)
),
points AS (
    SELECT id, st_setsrid(st_point(x, y), 32616) geom
      FROM coords
)
SELECT id, ST_Scale(st_buffer(geom, 1), ST_Point(4, 8), geom) geom
  FROM points
)


texts: Does this ArcGIS flow accumulation example have an error?

texts:
The lower right corner of the flow accumulation example in ArcGIS Desktop's documentation has a flow accumulation value of 2.

However, I only see one cell flowing into the lower right corner (the cell immediately above):

Is there some nuance to the algorithm that I'm missing that would result in the flow accumulation being equal to 2 at this cell?

I don't care about the software being used - just the example/algorithm.


texts: In a QGIS-project I have a shapefile consisting of polygons. However some of the polygons are double and thus laying over each other. Is there a way I can select all these 'conflicting' polygons so I can delete them from the shape?


texts: GDAL ogr2ogr OracleSpatial to Shape: Multipolygons with wholes/multiparts are missing

texts:
We're holding geodata in an Orcale database (11g) with spatial extension. I want to unload specific polygons and convert them into Shape for further local processing.

The table where I'm unloading from holds polygons/multipolygons: Oracle Spatial GeomTypes are 4003 and 4007 (4D-Polygon and 4D-Multipolygon - third and fourth axis are not used however).

For unloading and conversion I use this command within a python-script:

ogr2ogr.main(["", "-f", "ESRI Shapefile", result_path + query_result, connection_params,"-nlt", "MULTIPOLYGON", "-sql", the_sql_query_here])

Connection and convertion only seems to work for polygons with GeomType 4003 however:

Multipartpolygons or polygons with wholes (GeomType 4007) are completely ignored (white gaps). I already tried to rerun the command above with POLYGON or PROMOTE_TO_MULTI instead of MULTIPOLYGON. No change in the result though.

Any known caveats or workarounds regarding this problem?

My GDAL is v2.1.2

As requested in the comments I publish the geom-content of an object whose geometry couldn't be loaded:

MDSYS.SDO_GEOMETRY(4007,NULL,NULL,MDSYS.SDO_ELEM_INFO_ARRAY(1,1005,5,1,2,1,5,2,1,9,2,1,13,2,1,17,2,1,25,1005,5,25,2,1,29,2,1,33,2,1,37,2,1,41,2,1),MDSYS.SDO_ORDINATE_ARRAY(33433109.145,5711835.706,9999.9999,0,33433166.721,5711687.09,9999.9999,0,33433229.435,5711904.544,9999.9999,0,33433205.628,5711904.647,9999.9999,0,33433175.618,5711821.022,9999.9999,0,33433109.145,5711835.706,9999.9999,0,33433020.391,5711673.812,9999.9999,0,33433153.215,5711635.136,9999.9999,0,33433164.011,5711674.531,9999.9999,0,33433099.673,5711838.076,9999.9999,0,33433071.502,5711844.871,9999.9999,0,33433020.391,5711673.812,9999.9999,0))

Corresponding real life object is multipart and looks like this:

While the geometry doesn't get converted, the objects non-spatial attributes and values get fully converted to shape .


texts: Adding Altitude to a marker in WebGLEarth

texts:
Im trying to add an altitude to a marker in WebGLEarth but i cant figure out how to do it. I cannot find anything in documentation either. The idea is that my markers act as satellites so i need the altitude from the earth.

//this doesnt work
var satellite_marker1 = WE.marker([0, 0], {elevation: 260.0}).addTo(earth);

If a marker is the wrong thing to use for this i have no issue using something that will work.


texts: Passing parameters to PostGIS SQL queries in Python: weird behaviour


texts:
I use QGIS 2.14. I created a simple model using processing --> grafic modeler

In QGIS, I use python console and it works:

processing.runalg("modeler:model_name","input","output")

CanI running this model outside QGIS?


texts: How to create time series animation with unequal time interval in the data in ArcGIS? for eg I have population data for the years 2013,2012,2010,2004,1990,1950,1900, so how can I create time series animation with this kind of unequal time interval?


texts: OpenGeo Explorer for QGIS 2.6.1 - No module named ogr2ogrtopostgis

texts:
I'm trying to use the OpenGeo Explorer plugin for QGIS (v2.6.1) to import shapefiles to my PostGIS database. I can connect to the server, but when I try to import, I get the following error:
No module named ogr2ogrtopostgis
I've been unable to find any reference to this module, or discussion of this error.

Windows 7


texts: Adding shapefile as layer in ArcMap via Python scripting?

texts:
I am trying to automate various tasks in ArcGIS Desktop (using ArcMap generally) with Python, and I keep needing a way to add a shapefile to the current map. (And then do stuff to it, but that's another story).

The best I can do so far is to add a layer file to the current map, using the following ("addLayer" is a layer file object):

def AddLayerFromLayerFile(addLayer):
  import arcpy
  mxd = arcpy.mapping.MapDocument("CURRENT")
  df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
  arcpy.mapping.AddLayer(df, addLayer, "AUTO_ARRANGE")
  arcpy.RefreshActiveView()
  arcpy.RefreshTOC()
  del mxd, df, addLayer

However, my raw data is always going be shapefiles, so I need to be able to open them. (Equivantly: convert a shapefile to a layer file without opening it, but I'd prefer not to do that).


texts: In researching this problem, it appears that there is an issue with loading extensions in sqlite3 that requires a recompile of sqlite3 so SELECT load_extension("mod_spatialite.dylib") can be called.   I have not been able to succeed recompiling sqlite3 on the mac or found anyway to access a spatialite 4.3 database using Python 3.6 or greater on the MAC. Has anyone been successful with this?


texts: Writing to file geodatabase in R on Windows?

texts:
The title above basically states the question. I have looked around but haven't quite found anything yet. I see this driver:

https://gdal.org/drivers/vector/filegdb.html

I haven't seen any R implementation of it yet. I see that one has to add this plugin to GDAL: https://trac.osgeo.org/gdal/wiki/FileGDB

Undoubtedly I'm not the first person to want to do this. Is anyone aware of any attempts to implement this in R?


texts: Installation of GDAL 1.10 in Windows 8 and Ubuntu 12.04

texts: Is there an easy way to install GDAL/OGR 1.10 in Windows 8 and Ubuntu 12.04.


texts:
I'm adding some layers on a base map using OpenLayers 2. The aim is to create a layer with some features in it and later when user changes his prefered query, some other features need to be shown, so its the same layer but with different features. Therefore, I need to remove the layer and create it again with new features and then show it on the map.

I'm using the following code for defining and adding the layer:

var Selected_Buildings;
Selected_Buildings = new OpenLayers.Layer.Vector("MyLayer",  { eventListeners:{
            'featureselected':function(evt){
//some actions here
            },
            'featureunselected':function(evt){
//some other actions here
            }},
        strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "url to php code that runs query and return geo-data in kml format",
                format: new OpenLayers.Format.KML({
                    extractStyles: false,
                    extractAttributes: true,
                    maxDepth: 2})
            }),
                styleMap: new OpenLayers.StyleMap({
                    "default": new OpenLayers.Style(OpenLayers.Util.applyDefaults({
                       strokeWidth: 2,
                     strokeColor: "Gold",
                     fillColor: "Gray"
                    }, OpenLayers.Feature.Vector.style["default"])),
                    "select": new OpenLayers.Style({
                     fillColor: "Blue"
                    })
                })
            });

And the code for defining SelectControl:

var selectControl;
selectControl = new OpenLayers.Control.SelectFeature(
                [Selected_Buildings],
                {
                    clickout: true, toggle: false,
                    multiple: false, hover: false
                }
            );
                map.addControl(selectControl);
                selectControl.activate();

and then when user clicks on a button, I try the following code for removing my layer and then again the above code is executed to add new layer:

var mLayers = map.layers;
var remLayers = [];

for (var i=0; i<mLayers.length; i++) {
    if (mLayers[i].name == "MyLayer") {
        remLayers.push(mLayers[i]);
    }
}

for (var i2=0; i2<remLayers.length; i2++) {
    map.removeLayer(remLayers[i2]);
}

map.removeControl(selectControl);

I receive no error, but also the layer is not removed from the base map.

Any suggestion how to solve this problem?

In the above code MyLayer is the name of layer I would like to remove. selectControl is the name of the map control that was created before for the map.


texts: What does "Error: Layer not found: -1" from ArcGIS Silverlight API 2.4 mean?

texts:
We are using ESRI SL 2.4 API- and our application has one dynamic service and one feature layer. Only few of the users getting below error, not sure whether this is issue with Microsoft sliverlight version. Any quick help on this.

Error: Layer not found: -1

at ESRI.ArcGIS.Client.Layer.OnInitializationFailed(EventArgs e)
   at ESRI.ArcGIS.Client.Layer.Initialize()
   at ESRI.ArcGIS.Client.FeatureLayer.info_Error(Object sender, TaskFailedEventArgs e)
   at ESRI.ArcGIS.Client.FeatureService.FeatureLayerInfo.OnError(Exception err)
   at ESRI.ArcGIS.Client.FeatureService.FeatureLayerInfo.Execute_Completed(Object sender, DownloadStringCompletedEventArgs e)
   at System.Net.WebClient.OnDownloadStringCompleted(DownloadStringCompletedEventArgs e)
   at System.Net.WebClient.DownloadStringOperationCompleted(Object arg)


texts: V.net.connect result causing "not responding" error

texts:
I have a point shapefile layer (744 points) that I want to connect to a line shapefile layer (42 lines) using v.net.connect from GRASS 7. The process gives me the desired output, even if in wrong CRS (User:100000 instead of the CRS of the input layers) and as separate segments (1200+) instead of 42 single multi-line features as desired.
The main problem is that every operation I try to do with the newly created layer causes the software to go on "not responding" mode and eventually crash.

I tried multiple actions such as changing CRS, performing singlepart to multipart, adding a virtual field or even just opening the attribute table, but each of these operations causes the software to crash. (All other layers are still working perfectly fine).

Has anybody else experienced that or has any idea on what may be the problem and how to solve it?


texts: Is there an algorithm for transforming from CGCS2000 to WGS84?

texts:
Is there an algorithm that supports transforming from CGCS2000 coordinate to WGS84 coordinate? CGCS2000 is China Geodetic Coordinate System 2000.

Any advice would be thankful!


texts: ArcGIS & historical versions performance expectations

texts:
Short background: ArcGIS Server 10.1.1 Enterprise edition using MSSQL for our RDBMS.

Currently we utilize versioning only to support replication between 2 databases (edit to DEFAULT). Of late, I'm being asked more and more about the potential for producing 'snapshots' in time. Of course, ESRI supports historical based versions, but I get the sense that this is more for project-level management...not long term management, mostly because for performance reasons, you'll need to be reconciling versions, etc. We also conduct many many edits even on a given day; this DB is constantly updated.

Questions therefore are:

Am I interpreting the ESRI documentation correctly in terms of historical versions not being long-term data management strategies?
Are there other options I should be looking at to support this request for historical data modeling/tracking/archiving?
For those using ESRI's archiving: what has your performance been like? How long have you sustained historical versions for?

Thanks in advance!


texts: Calculating mean index value of raster for grid cell map using ArcMap

texts:
I have a raster index map and and also a polygon fishnet. I would like to assign mean index value into the grid cells.

How can I do that in ArcMap?


texts: Merge attributes using spatial join where Target and Join feature is the same

texts:
I have a single polyline feature class and i need to merge one of the attributes together based on a common identifier using a delimiter (i.e. |).

I was certain Spatial Join using the one feature class as both the Target and Join features would be the answer (I have tried various settings tried using ONE TO MANY, 'Identical To', join with pipe delimiter on 'RUFI_wlen').

But it doesn't produce the results. There are multiple rows per UFI in the results.

Using ESRI ArcMap 10.5.

This image is the input attribute table:

The desired result would be:

ROW 1:
1152, 03329_R M WILLIAMS WAY_RC__10|03492_AYR ST_RC__18|03580_AYR ST_RC__61
ROW 2:
1153, 03329_R M WILLIAMS WAY_RC__10| 03492_AYR ST_RC__18 |03580_AYR ST_RC__61
etc.


texts: Getting All Paired Shortest Paths in ArcGIS network analysis at once?

texts:
I have some difficulties on network analysis,I want to run shortest path algorithm,
 I  don't know how to use ArcGIS to reveal the shortest route between the appointed points on the map, and if it's possible, can I also use this software to record and calculate the number of the real road be taken as shortest route?

I have about 200 point to be the origins and destinations on the map, and I need the shortest route between these appointed points (200*200) and let them revealed on the map
And then,record and calculate the number of the real road be taken as shortest route.

but i can only have one shortest route between two points once a time if i use the toolbar, i am wondering how can i make all the shortest routes between 200 points one time? do i have to use any progamming method?


texts: Using OpenLayers with alternative map source to integrate Google Street View?

texts:
I am researching any possibilities of Using an alternative map source other than Google that I can integrate Google Street View.

Due to the drastic rate changes for Google Maps API I would like to use OpenLayers with Here Maps API and integrate only Street View from Google Maps API.

Is this possible?


texts: Does it make sense to create graticules using UTM meters?

texts:
Here's what I think I know, but I'm a geography newbie, so please correct any misconceptions I have, especially with regard to the term graticule.

A graticule is a specialized grid where the lines represent parallels and meridians.  If those lines are latitude and longitude, then this kind of graticule is called a "conjugate graticule".

This is what I think is the most common usage of the term "graticule".  However, I think it's just a common usage (and not a necessity) for the graticule's lines to be latitudes and longitudes.

My question is, if you have a UTM based map, isn't it possible (and more useful) to have a graticule in UTM coordinates (meters)?


texts:
I wanted to save a raster (1 band) file as a jpeg with:
xarray.rio.to_raster(save_path, compress='JPEG',
                         photometric='YCBCR' if xarray.sizes['band'] == 3 else None,
                         num_threads='all_cpus', tiled=True, dtype='uint8')

The nodata value originally is 55537.
I get the error:
error:  rasterio._err.CPLE_AppDefinedError: MissingRequired:TIFF directory is missing required "StripOffsets" field
I was wondering, why I get an Error for GTIFF while saving a JPEG, therefore I added the parameter:
driver='JPEG', which didn't change anything.
I already save the raster with:
    xarray.rio.to_raster(save_path, compress='zstd', zstd_level=1, num_threads='all_cpus', tiled=True, dtype='uint8', predictor=2)

and didn't encounter any error. What cold cause the error?


texts: Creating solid polyhedral surface from extruded 2D polygon using PostGIS

texts:
How can I create a solid polyhedral surface from extruded 2D polygon? I have created a closed polyhedral surface but ST_MakeSolid does not work for me.
Is there something wrong with the code below?
WITH cte AS (
    SELECT ST_MakeSolid(
        ST_Extrude(
            ST_Translate(
                ST_Force3D(
                    ST_GeomFromText('POLYGON ((-1 1,1 1,1 -1,-1 -1,-1 1))')
                ),
                0, 0, -1
            ),
            0, 0, 2
        )
    ) geom
)
SELECT
    ST_IsSolid(geom),
    ST_IsClosed(geom),
    ST_Volume(geom)
FROM cte

# OUTPUT  : false, true, 8
# EXPECTED: true, true, 8


texts:
Would anyone know the reason why this happens? and possibly a way to fix it. I am stumped

In the following code. The first addlayer works, but crashes on the second addlayer with a 'vshost32.exe has stopped working' error.

However when I comment out all the addlayers and make it to the for loop.  All three addlayers display correctly.

Why would the exact same code work in a for loop, but not when add a second, third and fourth layer individually.

Any ideas?

I have taken care of the error, however when I stop the program, I loose the data in Arcmap 10

// want to open shapefile
pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass("wa_body_" + sGeoname);
pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");

pFeatureLayer.FeatureClass = pFeatureClass;
pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
pFeatureLayer.Name = pDataset.Name;

pMap.AddLayer(pFeatureLayer);

pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass(QCForm.chkLstNSHNPolygons.Items[8].ToString()); // "wa_poly_" + sGeoname);
pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");
pFeatureLayer.FeatureClass = pFeatureClass;
pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
pFeatureLayer.Name = pDataset.Name;

pMap.AddLayer(pFeatureLayer);

pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass(QCForm.chkLstNSHNLines.Items[4].ToString());// "nshn_line_" + sGeoname);
pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");
pFeatureLayer.FeatureClass = pFeatureClass;
pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
pFeatureLayer.Name = pDataset.Name;

pMap.AddLayer(pFeatureLayer);

pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass(QCForm.chkLstNSHNPoints.Items[5].ToString()); // "nshn_point_" + sGeoname);
pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");
pFeatureLayer.FeatureClass = pFeatureClass;
pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
pFeatureLayer.Name = pDataset.Name;

pMap.AddLayer(pFeatureLayer);

pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass(QCForm.chkLstNSHNJunctions.Items[4].ToString()); // "nshn_junc_" + sGeoname);
pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");
 pFeatureLayer.FeatureClass = pFeatureClass;
 pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
 pFeatureLayer.Name = pDataset.Name;

 pMap.AddLayer(pFeatureLayer);

 //Add topology layers
 i = 0;

 for (i = 0; i <= lstFeatures.Items.Count - 1; i++)
 {
     pFeatureClass = pLocalFeatureWorkspace.OpenFeatureClass(lstFeatures.Items[i].ToString());
     pFeatureLayer = (ESRI.ArcGIS.Carto.IFeatureLayer)pObjectFactory.Create("esriCarto.FeatureLayer");
     pFeatureLayer.FeatureClass = pFeatureClass;
     pDataset = (ESRI.ArcGIS.Geodatabase.IDataset)pFeatureClass;
     pFeatureLayer.Name = pDataset.Name;

     pMap.AddLayer(pFeatureLayer);
 }

}


texts:
I've been trying to import .cel files as symbology in QGIS. I was recommended to change those .cel files into .dxf, which I did. However, in order for the conversion from dxf to xml to work, I need the XML schema sample for QGIS symbology. I've exported existing QGIS symbology as an xml file. However, editing those exports to include the symbols I want to use, doesn't work.
What should I do?


texts: Route starting from closest node along edge

texts:
I am trying to find a way to route from the closest point on a linestring to my current location (lat, long). So far I am able to get the shortest path but it starts from the very beginning of the linestring (aka source). I am using prg_trsp
http://docs.pgrouting.org/2.0/en/src/trsp/doc/index.html
because it has a feature to specify the starting position along the linestring. I am able to correctly calculate the distance along the linestring and pass the values to the function but cannot figure out how to use the results from the function (pgr_costResult[]) to specify where the route should start (partially along the closest linestring).

I have a feeling I am doing something wrong with the join when I go to join the results from the routing algorithm to my edge table to get the geometry because when I join it uses the edge table's full geometry and not segments. Although, looking at the documentation, I don't see where you get a returned segment from the routing function.

Below is a screenshot of what I am trying to do (red line) and what I have (blue line) the point is the current location. The red line comes from using the pgrouting plugin in qgis with the trsp(edge) selection.

SELECT st_linemerge(edgeTable.geom_way) FROM pgr_trsp('SELECT id, source, target, cost FROM edgeTable',
    (SELECT id FROM origin),
    (SELECT * FROM sourcePos),
    (SELECT id FROM destination),
    (SELECT * FROM destPos),
    false, false) AS shortestPath
    JOIN edgeTable ON shortestPath.id2 = edgeTable.id;

origin is the id of the starting route

sourcePos is how far along the linestring to offset

destination is the id of the end linestring

destPos is the fraction of the end linestring

all as specified here: http://docs.pgrouting.org/2.0/en/src/trsp/doc/index.html


texts:
I am trying to set privileges in my database for users so they can only UPDATE, INSERT, DELETE and, of course, SELECT on geographic tables.

Since I am using some SERIAL types for my IDs, I would prefer to not even give the users the possibility to manually edit this field. So I am defining column-wise privileges. It works fine on the UPDATE part, but the button in QGIS allowing the creation of a new object stays greyed. The only way it seems to be working is when I do not specify any field in the privileges definition. Even when I choose all the fields it does not work (even though I would have thought that not specifying any column and specifying all of them would be the same).

It seems to be something that I don't understand, or there is a limitation in the privilege definition for proper interaction between the DB and QGIS. Does anyone have any piece of information or advice which could help me understand what is happening, and/or (even better) help me achieve my goal?

I can always deal with that by setting the field as non-editable in the style definition, but since anyone can set it as they please, I would prefer a more secure alternative.

Running QGIS 2.14, PostGIS 2.3 for PostgreSQL 9.5.


texts: I am trying to make web map with the help of Qgis2leaf plug in. i added all layer to be plot but when i clicked on OK after setting all parameters i am getting Python coding error. Here i am attaching my error snapshot. Is it related to QGIS version? or related to python plugin for QGIS2Leaf.


texts: PlanetScope Image Delivery to Google Earth Engine

texts: I am trying to order PlanetScope data from Planet Explorer and have it delivered directly as an asset to my image collection in Google Earth Engine, but it is unclear how to do this.  I have read and implemented everything the documentation and watched the one youtube video on this subject, yet I still cannot figure out how to send the order (of my selected scenes) to my account.  This is my first time working with Planet imagery, so perhaps this task is more straightforward than what I am experiencing.


texts: gdalwarp for clipping images using multiple polygon in a bash shell for loop

texts:
C:\OSGeo4W64\bin\gdalwarp.exe -dstnodata 0 -q -cutline D:\hillsborough_Images_work\split_Polygon\120-140k\Hillborough_FINAL_140000_region_IGD_120570003d29ed0dc.shp -crop_to_cutline -dstalpha -tr 0.500000000000000 0.500000000000000 -co "COMPRESS=LZW" D:\hillsborough_Images_work\HIL2017\mrsid.vrt D:\hillsborough_Images_work\Clipped_Images\120570003d29ed0dc.png

C:\OSGeo4W64\bin\gdalwarp.exe -dstnodata 0 -q -cutline D:\hillsborough_Images_work\split_Polygon\120-140k\Hillborough_FINAL_140000_region_IGD_120570004819de834.shp -crop_to_cutline -dstalpha -tr 0.500000000000000 0.500000000000000 -co "COMPRESS=LZW" D:\hillsborough_Images_work\HIL2017\mrsid.vrt D:\hillsborough_Images_work\Clipped_Images\120570004819de834.png

C:\OSGeo4W64\bin\gdalwarp.exe -dstnodata 0 -q -cutline D:\hillsborough_Images_work\split_Polygon\120-140k\Hillborough_FINAL_140000_region_IGD_120570006c12a2d77.shp -crop_to_cutline -dstalpha -tr 0.500000000000000 0.500000000000000 -co "COMPRESS=LZW" D:\hillsborough_Images_work\HIL2017\mrsid.vrt D:\hillsborough_Images_work\Clipped_Images\120570006c12a2d77.png

C:\OSGeo4W64\bin\gdalwarp.exe -dstnodata 0 -q -cutline D:\hillsborough_Images_work\split_Polygon\120-140k\Hillborough_FINAL_140000_region_IGD_120570008cca2ea3a.shp -crop_to_cutline -dstalpha -tr 0.500000000000000 0.500000000000000 -co "COMPRESS=LZW" D:\hillsborough_Images_work\HIL2017\mrsid.vrt D:\hillsborough_Images_work\Clipped_Images\120570008cca2ea3a.png

This worked fine but I wanted to do this task using for loop and I writed script for this but I'm doing something wrong. Can you please help me to do it correct?

#!/bin/sh
#    "shp" - folder for shapefiles
#    "outputraster" - folder for output rasters
cd /home/user/cliprasters/
#  delete previous files
rm D:\hillsborough_Images_work\temp\*.png
#  run
for fname1 in D:\hillsborough_Images_work\split_Polygon\0-20k\*.shp do
  fname2="$(basename $fname1 .shp)"
  C:\OSGeo4W64\bin\gdalwarp.exe -dstnodata 0 -q -cutline D:\hillsborough_Images_work\split_Polygon\0-20k\$fname2.shp -crop_to_cutline -dstalpha -tr 0.500000000000000 0.500000000000000 -co "COMPRESS=LZW" D:\hillsborough_Images_work\HIL2017\mrsid.vrt D:\hillsborough_Images_work\temp\$fname2.png
done


texts: QgsFeature.setAttributes(attr) blocks access to changing attribute value

texts:
I have a problem with QgsFeature.setAttributes() method, I have tried something like the code below:

some_layer = iface.activeLayer()

fields = some_layer.fields()

extra_fields = [QgsField('field1',QVariant.String), QgsField('field2', QVariant.String)]

for field in extra_fields:
    fields.append(field)

some_fts = some_layer.getFeatures()

for some_ft in some_fts:
    attrs = some_ft.attributes()
    new_feature = QgsFeature(fields)
    new_feature.setAttributes(attrs)
    new_feature['field1'] = 'some text'

When this is done, the new feature 'field1' attribute has value NULL.
When I out-comment the line:

#new_feature.setAttributes(attrs)

the new feature 'field1' attribute has value 'some text', ofcourse in that case all the other attributes are nulls.

Does anyone know why accessing the attributes, after setAttributes method has been used, doesn't work?


texts:
Is there any way so i can add symbols like ( <, > ,- ) in QGIS polygons/vectors so I can define for a example in which way is "this pipe of drainage leaking water"
Working in QGIS 3.2 Bonn


texts: How to extract features from a vector layer?

texts:
I have a vector file containing the river lines of an area. I want to extract a specific water line from the vector. I used:
v.db.select but it only displayed:

   cat|F_CODE_DES|HYC_DESCRI|NAM|ISO|NAME_0
    1|River/Stream|Perennial/Permanent|UNK|PHL|Philippines
    2|River/Stream|Perennial/Permanent|UNK|PHL|Philippines
    3|River/Stream|Perennial/Permanent|UNK|PHL|Philippines
...

then i used v.extract to select a specific area and it displayed this error:

v.extract input=vector_waterLines@PERMANENT output=vector_1 where=1
DBMI-DBF driver error:
ERROR: Unable select records from table <vector_waterLines>


texts: How can I scale LiDAR intensity values to 8-bit, in order to apply the Wing2015 snag algorithm using segment_snags in lidR?

texts: I'm trying to identify snags using LiDAR stored in las format. The Wing 2015 algorithm used in segment_snags requires the values be in the range from 0 to 255. Is it possible to edit and scale intensity values (not by adding a new attribute)? Alternatively, is it possible to edit the segment_snags tool in lidR?


texts:
I am trying to make a world population density Choropleth map (country level).

The data I am using is publicly available on World Development Indicators website.

Pretty much, I have an excel file with one column for country name and the other column for population density.

I want to create a simple world map but I have very minimal experience in using QGIS (I have a mac). I tried using google maps but rather than shading, it places markers on the maps.

I want something simple that looks like this: http://1.bp.blogspot.com/-Gl7YYBhOuwI/TmQPm7qrnuI/AAAAAAAAAAU/lslkxp8OZ4w/s1600/choropleth+map+world.jpg

and has a nice scale which shows what the max and min population density are.

Can you guys direct me to any good resources on this? (Keep in mind that I have very little experience with QGIS)

What would be the best way to make such a map?

P.S. I know I am a noob :( . Can someone please direct me towards some good introductory resources for using QGIS/ArcGIS and mapping software in general? This will really help me with my job.


texts: Jagged lines when increasing line width in QGIS

texts:
I have a curved multipart line that when the line width is increased the edges become jagged. I'd like the line to be smooth if possible. I am using QGIS 2.19 Las Palmas. Is there any wat I can do this?

Apologies for my lack of correct terminology - I am new to QGIS.


texts: How to use the R function: raster::extract with polygon which partly covers raster?

texts:
I have a raster file which is only partly covered by a shapefile with polygons.

Now I want to apply raster::extract in order to get the mean value for each polygon based on the rasterfile - But some features of the shapefile are complete outside of the rasterfile.

If I'm using snow clusters in order to speed the extraction up - I get the following error:
R raster::extract Error in .local(x, y, ...) : cluster error at polygon: 1

Do you have an idea, how to resolve this issue?


texts: Applying equation to a Numpy array while preserving TIFF metadata (coordinates)

texts:
I wish to load a .tif image into a numpy array, apply a formula and then spit the image back out with the same geo/nodata metadata information.

So far I have created this, which does the processing but loses all metadata in the process and also applies the formula to nodata values:

from PIL import Image
import numpy as np

img = Image.open('ndvi.tif') #import image array
ndviArray = np.array(img) #convert to numpy array

def ndvi2bio(x): #function to convert from ndvi to BIO
    return 400*(x)+680

vfunc = np.vectorize(ndvi2bio) #creating function function array...?

bioArray = vfunc(ndviArray) #applying the bio array to the formula array to give results.

#turn back into image...
bioImg = Image.fromarray(bioArray)

bioImg.save('BIOout.tif')

I've been reading into rasterio, I think I can replace the PIL image module with rasterio which will be more appropriate in this situation and allow me to keep my metadata from the TIFF.

How can I achieve this? So far I am struggling to find a tutorial which does this, it seems most others use several other modules. Can it be done with just rasterio and numpy? (img(+meta) -> numpy array -> img(+meta from original data))


texts:
Possible Duplicate:
  How to skip drag event to stop a feature being dragged

I've created a map and added markers to baselayer on click by using this code :

vectorLayer = new OpenLayers.Layer.Vector("Vectors");
        vectorLayer.id = "Vectors";
        map.addLayer(vectorLayer);

        map.events.register("click", map, function (e) {

            //gets pixels on mouse position
            position = map.getLonLatFromPixel(e.xy);
            latitude = position.lat;
            longitude =position.lon;

            //Adding markers to vectorLayer
            marker = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Point(longitude, latitude),
                { draggable: true },
                { externalGraphic: 'http://www.openlayers.org/dev/img/marker-blue.png',
                graphicHeight: 25, graphicWidth: 21,
                graphicXOffset: -((size.w / 2) - 2), graphicYOffset: -(size.h - 2) });

                // function markersIdGenerator() will return a random id for the marker
                marker.id = markersIdGenerator();
                // Adding marker to vector layer
                vectorLayer.addFeatures(marker);
            });

I've added DragFeature control to drag markers is there any way to fix a marker ?

I've tried like this:

draggable=true;
    dragmarker = new OpenLayers.Control.DragFeature(vectorLayer, {
        autoActivate: true,
        onEnter: function (marker) {
            if (dragChecker(marker) == true) {
                draggable = true;
                dragmarker.activate();
            }
            else {
                dragmarker.deactivate();
                draggable = false;
                return;
            }
        },
        onLeave: function (marker) {
            if (draggable == false)
            { return; }
            else {
                Log("draggable value => "); Log(draggable);
                dragmarker.activate();
            }
        },
        onComplete: function (marker, pixel) {
            if (draggable == false)
            { return; }
            else {
                var fixer = confirm("Are you want to Fix this marker");
                if (fixer == true) {
                    marker.data.draggable = false;
                    console.log("Fixed");
                }
                else {
                    marker.data.draggable = true;
                    console.log("Not Fixed");
                }
            }
        }
    });
    map.addControl(dragmarker);

    //  this function will return whether the marker is dragable or not
    //  by checking marker.data.draggable
    function dragChecker(marker)
    {
        if (marker.data.draggable == true)
        { return true; }
        else
        { return false; }
    }


texts: How to use expression to modify feature class name in ModelBuilder output?

texts:
I am using ModeBuilder to project many raster files from one coordinate system to another.  I'd like the output file names to be different than the input.

The file names currently end in "utm83" but I want to remove that and add "NAD83".  I know how to use the %name% variable to keep the same name, but I'm not sure how to alter that for my desired effect.

It seems like I should be able to use the "Calculate Value" with an expression like LEFT(%name%, -5) to get this result but I can't figure it out.

Here is what my model looks like:


texts:
I have a SRTM DEM and a DEM derived from satellite (Cartosat). Both the DEMs are registered perfectly. Now I want to carry out a DEM 'fusion' process to derive a higher quality output DEM.

I am not using the fusion software. The two dems are not having any holes. I need tools to fuse two DEMS to generate a better quality DEM. (Weighted averaging, sparse representation or frequency based filtering )

I want to carry out in arc gis.


texts: How to ensure different colors for surrounding polygons in ArcGIS?

texts: I can symbolize my polygon in ArcGIS by categories. I tried for every field. But unfortunately in few polygon it showing same color those are just beside of same color polygon. I have to categorized 572 polygon. So it will be a hectic task if I have to do this manually. Is there any way to categorized polygons which are not identical in color with the polygons its surrounded?


texts: How to create/select a unique Identifier OID in case I have a spatial view and none of the fields or their combination can be unique

texts: How to create/select a unique Identifier OID in case I have a table view in oracle and none of the fields or their combination can be unique. The process of selecting OID is described here and here.


texts: Take centroid of multipolygons and use it as center of a circle

texts:
I'm working with QGIS (or rather discovering QGIS, should I say) and I have a data structure that looks like this:

PLACE_NAME      ST_ASTEXT(geom)
     x        MULTIPOLYGON1(...)
     x        MULTIPOLYGON2(...)
     y        MULTIPOLYGON3(...)
     y        MULTIPOLYGON4(...)
     y        MULTIPOLYGON5(...)
     x        MULTIPOLYGON6(...)
    ...              ...

For a certain place (say "x"), I would like to compute the centroid of all its multipolygons taken as a whole. Then, I need to use this centroid to find all the multipolygons that are located within a circle whose center is the above-mentioned centroid and whose radius is 10 kilometers. In a mix of pseudo SQL and personal code, this looks like this:

SELECT centroid(average(geom)) as UComP -- Unique Centroid of my Polygons
FROM my_table
WHERE PLACE_NAME = 'x'

And then:

SELECT *
FROM my_table
WHERE ST_Within(geom, ST_Transform(ST_Buffer(ST_Transform(ST_SetSRID(ST_MakePoint(ST_X(UComP), ST_Y(UComP)), 4326), 3857), 10000), 4326)) = true)

Unfortunately, I'm facing some problems with the beginning, i.e. the calculation of the centroid. It seems that PostGIS doesn't really like the average function called on a geometry.

In addition to a solution to the above problem, I'd ideally need a unique query instead of two distinct ones.


texts: can't calculate numeric field with ogr python

texts:
I have created two functions: one to add a field to a layer and one to calculate a field. I am having issues with the calculate field function when trying to calculate values for a numeric column. I do not get an error, but the column is blank after the code runs (it was blank to begin with). I've tried both OFTInteger and OFTReal as field types with the same result. Here are the functions:

def add_field(lyr, field_name, field_type):
    field = ogr.FieldDefn(field_name, field_type)
    field.SetWidth(50)
    lyr.CreateField(field)

def calc_field(calc_lyr, field, exp):
    for feature in calc_lyr:
        geom = feature.GetGeometryRef()
        area = geom.GetArea()
        feature.SetField(field, exp)
        calc_lyr.SetFeature(feature)

and here are the function calls to add and calculate a numeric field:

add_field(dxf_lyr, 'color', ogr.OFTInteger)
calc_field(dxf_lyr, 'color', 150)

The field gets added, but the values don't get caculated. What am I missing?


texts: Creating Shapefile/ Feature class from Grid CSV using ArcGIS for Desktop?

texts:
I have csv file which contains data in grid form as:

I have created a python Script that convert it into below form:

Is there any tool available within ArcGIS for Desktop that can i used to do my job rather than creating a new script/tool and create vector points geometry in those feature classes.


texts: Vertical Datum transformation using Pyproj

texts:
I have been using Pyproj for some time now doing only simple projection and transformation. I am totally new to vertical datum and differences in GEOID among Coordinate systems.
I have been trying to do vertical transformation using pyproj but I feel like I am missing something.
For example I tried this :

import pyproj
lat = 43.70012234
lng = -79.41629234
z = 100
wgs84 = pyproj.Proj("+init=EPSG:4326")
NAD83NDV88 = pyproj.Proj("+init=EPSG:5498")
results = pyproj.transform(wgs84, NAD83NDV88, lng, lat, z)

which returns : (43.70012234, -79.41629234, 100.0).

My logic was to take a WGS84 coordinate with an ellipsoid elevation and try to convert that into NAD83 NADV88 (epsg:5498). But it seems like it is not that simple.

I found some EPSG codes for transformation related to this question like : EPSG:3858 +> WGS 84 to EGM2008 height (1) but I haven't been able to figure out how to use those codes in Pyproj. I am missing some pieces of the puzzle!

So to summarize can this be done?


texts: What's the best way to visualize hundreds of GPX files?

texts:
I'm an amateur GIS enthusiast who maps my runs, hikes and urban explorations all over the place. I have accumulated hundreds of GPX files which I enjoy visualizing in different programs - but all of them fall short for one reason or another.

Googling just brings me to a lot of nonsense about interoperability between different run tracking programs and google maps.

Do any of you experts have a suggestion for a tool that will let me visualize all my GPX tracks in one world map? The closest thing I came across was umap but it had some annoying artifacts with certain files.


texts:
When I tried to create a buffer layer around a polyline layer. But I got the error as follows:

The input was outside a +-360 longitude range or outside +-90 latitude.

The input is : dev_db.DBO.Segment

The output is : D:\local\SegBuffer.shp

The buffer size : 10 meters

The error occurs regardless of the buffer size. Any hint on it?

I am using arcmap 10.2.2.


texts: Non-Research Transportation Raster applications

texts:
My experience with rasters is pretty basic. 99.95% of what I do involves vector and transportation modelling. In research, there are numerous transportation applications using rasters and imagery however, these applications although realistic and great, seem to me as a bad representation of the state of practice of Transportation and rasters/imagery.

What is the state of the practice on using imagery and raster analysis in transportation?


texts: How do I cluster XY point data? (ArcGIS 9.3)

texts:
I have a whole bunch of icons displaying on a map, I want these clustered together as was suggested to me for a different question here...  ArcGIS 9.3. How do I increase Picture Marker Symbol performance?

A link to an example was provided but I'm not really sure where to begin.  I'll preface this by saying I'm a programmer way before I'm a gis person.

What I want to know, in total noob terms, is how can I cluster my icons for XY events on the map?

Thanks!

EDIT: Forgot to mention, this is for a C# web ADF GIS portal... If it makes any difference.

Is this even possible?  I can't seem to find anything for doing this with 9.3 adf.


texts: Fill raster with closest values from other raster with rastercalculator

texts:
I have two raster:

R1 --> Extension raster ("RB_100" in the picture). This raster is used as expected extension for data and contains just 1 where I need data, and 0 where I don't.

R2 --> Precipitation raster ("Precipitation" in the picture). With values from 1 to 10.

Both rasters are already in the same extension and resolution. What I need is to generate a raster with data of R2 in every pixel where R1 is equal to 1 (grey pixels).

In other words, I need to create a raster where those grey pixels contain the closest value from R2 ("Precipitacion").


texts: Computing the harmonic mean of an image collection in the Google Earth Engine

texts: How can I compute the harmonic instead of arithmetic mean for a series of images in Google Earth Engine? There is a build-in function to compute the arithmetic mean and median, but there is no harmonic/geometric mean available. I tried to build my own function by making use of the reduce function, but with this I did not succeed. Any ideas/examples on how to proceed?


texts: QGIS Crashes when I connect to QGIS Cloud

texts:
QGIS crashes when I try to connect to QGIS cloud. I just signed up for the cloud service - still testing out the free version.

When I try to connect through the plug-in panel (I use no proxy server) the "wheel" on my mac starts spinning, QGIS stops responding, and I have to force-quit it.

I have tried it with and without the regular Mac firewall, makes no difference.

You can close this it's solved. When you have a problem like this you don't know what's causing it. So it wasn't off-topic. Anyway I solved it myself by getting a new router.


texts: How is distance between zip codes calculated?

texts: Specifically, where can I find a database of zip codes, their lat/lon centers, and intervening "obstacles" (for example, zip codes that are geographically close but separated by the Chesapeake Bay).


texts: Kriging with barriers in ArcGIS Desktop?

texts:
I am trying to do interpolation in a coastal area with data points in several rivers emptying into a bay and on opposite sides of peninsulas.

I have created a coastline layer to use it as a barrier but cannot figure out how to do this in ArcGIS 10.

I do not see a barriers object anywhere in any of the kriging tools.

I tried just doing the kriging then clipping it to the land, but I get unrealistic results where data points close together but separated by land create problems.

I see online documentation for doing this in geoprocessing scripts in 9.3 but nothing similar for 10.


texts: Make FeatureCollection from Features in Google Earth Engine

texts:
I try to add a list of Features into a FeatureCollection with ee.FeatureCollection but it raises the following Error: *

"FeatureCollection (Error)
Collection, argument 'features': Invalid type. Expected: List<Element>. Actual: List<List<Object>>"

I don't understand the difference between Object and Element. It is necessary to convert the single Features Object into an single Feature Element?

var s2aprop = ee.FeatureCollection([
ee.Feature(s2a_p.aggregate_array("CLOUDY_PIXEL_PERCENTAGE"),{name:'Cloudy_Percentage'}),
ee.Feature(s2a_p.aggregate_array("PRODUCT_ID"),{name:'Product_Id'}),
ee.Feature(s2a_p.aggregate_array("SENSING_ORBIT_DIRECTION"),{name:'Orbit_Direction'}),
ee.Feature(s2a_p.aggregate_array("SENSING_ORBIT_NUMBER"),{name:'Orbit_Number'}),
ee.Feature(s2a_p.aggregate_array("SNOW_ICE_PERCENTAGE"),{name:'Snow_Ice_Percentage'}),
ee.Feature(s2a_p.aggregate_array("VEGETATION_PERCENTAGE"),{name:'Vegetation_Percentage'}),
ee.Feature(s2a_p.aggregate_array("WATER_PERCENTAGE"),{name:'Water_Percentage'}),
ee.Feature(s2a_p.aggregate_array("UNCLASSIFIED_PERCENTAGE"),{name:'Uncalssified_Percentage'}),
ee.Feature(s2a_p.aggregate_array("SPACECRAFT_NAME"),{name:'Spacecraft'})
])

var s2aprops = ee.FeatureCollection(s2aprop)
print (s2aprop)

https://code.earthengine.google.com/ca44d644303c7212c4a11e4c00f752a9
Line 161-174


texts:
I am trying a convert a raster dataset [a simple 4 x 4] into a numpy array using RasterToNumPyArray. Here is my sample raster dataset:

RasterToNumPyArray however, returns zero for Null values:
inRaster = "Slope_NULL"

# Create a Raster Object and get its lowerleft
aRObject = arcpy.Raster(inRaster)
lowerLeft = aRObject.extent.lowerLeft

print aRObject.noDataValue

# Create a numpy array of the Raster Object
anArray = arcpy.RasterToNumPyArray(aRObject)
print anArray

and I get:

I am wondering if there is a way that I could get None or NaN instead of zero in my array?
I know I can change the Nodata to a value using nodata_to_value, however, RasterToNumPyArray syntax is explaining that the default value is None, and it is important for me to have None and not zeros in my array.
How can I have None or NaN [I am not sure which one is better supported in numpy] in my array?


texts:
Im trying to convert shapefile from one CRS to another in QGIS. On the layer menu, I right click on the shapefile, save it as with other name and there I change the CRS to the desired one (WGS 84). Then I check the properties of the file and under the tab of METADATA and it appers under Layer Spatial Reference System the desired CRS but when I project it, its not on the desired CRS. I also can see on the lower central coordinate display on the main screen that its not in WGS84 (because of the values). Any idea whats happening?

Thanks.


texts: The red box that refers to other dataframes will only show part of the box (only show over land but not water, see image)

texts: the red box refers to other dataframes but it's only showing over land and not over the ocea, ie, the box is an incomplete square


texts: How to setup PostGIS for use with multiple machines?

texts:
I'm a beginner user of postgresql/postgis. I am starting to use Linux more and more but still rely on windows (on different machines). How would I setup a postgis database that I could work with under different systems?
Thanks,
Eelke


texts:
I run a portable Linux box with an open street map tile server to provide maps for amateur radio emergency communications using APRS.

If needed to, I want a way to convert an address to a latitude and longitude where I do not have Internet connectvity.

Are there any Linux (in particular Ubuntu) geocoder utility's which will convert an address into lat/long?


texts: Transferring one field name from text file to new field name in feature class?

texts:
I am working on a script that  involves incidents of fire. I'm taking my values from a text file that is comma delimited. The first line of the text file contains the field names (Latitude,Longitude,Confidence). My goal is to create a new feature class that will take the last field name from my text file (e.g. Confidence) and use this as the field name in my new feature class. I have succeeded in getting arcgis to read the last field name from my text file so that a new feature class is created with the field names: "ID, SHAPE, Confidence", however, my feature class is empty- my for loops seems to stop after reading the first line. I believe that I have the right components, I'm just not sure what order to put them in.
The part that confuses me the most is, I start reading the text file in a for loop (which I think is pretty normal) and my instinct is to split this text file and, using indexing, point to the last field name "Confidence" using a variable. Once I do this, I'm not sure how to make my insert cursor and Add Field tool (which are outside of the for loop) recognize this variable that points to "Confidence". Any suggestions?

Below is my code:

import arcpy
work=r"C:\Users\Cara\Documents\WildlandFires.mdb"
arcpy.env.workspace=work

arcpy.env.overwriteOutput=True

iFile= r"C:\Users\Cara\Documents\NorthAmericaWildfires_2007275.txt"
with open(iFile, 'r') as f:
    lstFires = f.readlines()

#new feature class name
newFC= fireIncidents
fc=newFC

#user input mininum confidence threshold
confidenceThreshold=99
minThresh=int(confidenceThreshold)

lstFires = f.readlines()

#Create New Feature Class

arcpy.CreateFeatureclass_management(work, fc, "POINT", "", "", "", "", "",     "", "", "")

#Loop through file line by line and populate new rows in feature class
#Create Insert Cursor to use for inserting new rows into feature class

fields = ["SHAPE@", fieldName]
cur = arcpy.da.InsertCursor (fc, fields)

#add field
arcpy.AddField_management(fc, fieldName , "SHORT")
print "created fields"

cntr = 0 #counter variable used to print out progress of script

for fire in lstFires:

    if 'Latitude' in fire: #read line, use split function
        fn = fire.split(",")
        fieldName = fn[2]

        continue

    lstValues =  fire.split(",")#split out values into separate values

    if int(lstValues[2]) > minThresh: # if confidence value is greater than     user defined confidenc value, continue with code

        #create point object that will hold x,y coordinate of fire
        pnt = arcpy.CreateObject("Point")

        #pull out individual values using list indexing
        latitude = float(lstValues[0])
        longitude = float(lstValues[1])
        confid = int(lstValues[2])
        print "I can pull out values"
        cntr=cntr + 1

        #assign x, y properties of Point object to long and lat values
        pnt.X = longitude
        pnt.Y = latitude

        #Create new row for feature class
        #assigns new Point object to "SHAPE" field and confidence value to     "CONFIDENCEVALUE" field
        row = [pnt, confid]
        print "row's OK"
        #Insert new row to feature class
        cur.insertRow(row)
        print "Record # " + str(cntr) + " written to feature class"

    #release cursor lock on feature class
del cur


texts: Using custom behaviour on right click and actionPan on left click in QGIS plugins?

texts:
I implemented a maptool that draws a red cross on the point that is being right-clicked on. There are good how-tos, such as this topic : Programatically check for mouse click in PyQGIS?

But I happen to completely lose the left-click behaviour. This means that I cannot navigate through my layer with left click. I used a little trick ("iface.actionPan().trigger()") to get the pan action anyway, but it sets the mapTool to another object, an instance of mapToolPan. This later on causes various bugs, because I am calling methods that belong to my homemade PointTool class on a mapToolPan instance.

If the user right clicks on QGIS, the toolPan is set to my homemade PointTool class again.

Here is the code where I need help (I only implemented the canvasReleaseEvent because I only use this event) :

class PointTool(QgsMapTool):

def __init__(self, canvas):
    QgsMapTool.__init__(self, canvas)
    self.canvas = canvas

def canvasReleaseEvent(self, event):
    """ Checks which click was used : if left, triggers the pan action, else, adds point. """

    if event.button() == QtCore.Qt.RightButton:
        # The right click triggers the point selection
        x = event.pos().x()
        y = event.pos().y()
        self.clickedPoint = self.canvas.getCoordinateTransform().toMapCoordinates(x, y)
        self.showPoint(self.clickedPoint)
        self.canvas.refresh()

    elif event.button() == QtCore.Qt.LeftButton:
        print("release clicked left")
        # The left click triggers the pan action
        iface.actionPan().trigger()

def showPoint(self, point):
    """ Adds point in qgis main window """

    # Using a vertex marker to show the selected point as a red cross
    m = QgsVertexMarker(self.canvas)
    m.setCenter(point)
    m.setColor(QColor(255,0,0))
    m.setIconSize(10)
    m.setIconType(QgsVertexMarker.ICON_X)

    self.canvas.refresh()

How can I keep my right-click behaviour AND still have normal left-click behaviour?


texts: Default attribute when merging two polygons in ArcMap?

texts: When merging two polygons in ArcMap 10.3, it appears that the new polygon always inherits the attributes of the larger (by area) of the two input polygons, unless otherwise specified. Is this by design, or by coincidence?


texts: Exporting spatial object in mysql to shapefile with spatial reference

texts:
Info: Win7, MySQL 5.x and OSGEO4W 1.11

To import a shapeafile to MySQL I use this ogr2ogr command:

ogr2ogr -f MySQL MySQL:gis,user=root,password=password C:\posto_pianco.shp -nln pianco_post -a_srs EPSG:29194 -update -overwrite -lco engine=MYISAM

To export from MySQL to shapefile I use this:

ogr2ogr out_shapes MYSQL:gis,user=root,password=password

The Problem: The shapefile files exported have no spatial reference, how to resolve this?

Summary of layer:

Layer name: pianco_post
Geometry: Point
Feature Count: 1
Extent: (-37.928100, -7.204400) - (-37.928100, -7.204400)
Layer SRS WKT:
PROJCS["SAD69 / UTM zone 24S",
    GEOGCS["SAD69",
        DATUM["South_American_Datum_1969",
            SPHEROID["GRS 1967 Modified",6378160,298.25,
                AUTHORITY["EPSG","7050"]],
            TOWGS84[-57,1,-41,0,0,0,0],
            AUTHORITY["EPSG","6618"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4618"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","29194"]]
FID Column = OGR_FID
Geometry Column = SHAPE
objectid: Real (9.0)
codigo: Real (9.0)
nome: String (35.0)
altitude: Real (19.11)
longitude: Real (19.11)
latitude: Real (19.11)
disponibil: String (31.0)
area_drena: Real (9.0)


texts:
I know that it is possible to link databases from Access to QGIS via the ODBC Data Source administrator.

I am curious if it would be possible to create a python code/plugin that would allow the user to do the following:

Search for a specific piece of data
Create a buffer of x amount of distance from that one point
Select the data that intersect with that buffer

( I am sure the first 3 steps in this are possible)

Take the selected data and use it to spit out a report in Microsoft access via the ODBC Data Source Administrator.

Any thoughts?