texts: How To Order Labels from Several Layers in Mapserver?

texts:
I have an Mapserver project with several layers. Some of the layers have labels.

My problem is that labels from one layer overlap the labels from another layer. I want the labels to change vertical order (Label 2 should be rendered over Label 1 instead of Label 1 over Label 2).

I've tried to play around with "Force True/False" and "Priority," but either the Label 1 is above Label 2 or Label 1 isn't displayed at all.

The issue is that Label 1 has to be displayed certainly on the given x/y-point while Label 2 has to be dynamically positioned, because it is labeling a polygon and should be displayed in every zoom resolution and every image section.

The polygon which is described by Label 1 is overlayed by the Polygon of Label 2 (but view-able due to Opacity Settings). Label 1 itself isn't overlayed or overlapped by anything.

I hope I've described my problem clearly - does somebody know how to do that?


texts: How can I always run an algorithm with standard input shapefiles from mapcanvas in QGIS 2.0.1?

texts:
I created a model in the graphical modeler and I want to run that model always with two specific polygon shapefiles from my project, e.g. pst and astenot. Whenever I run my model, it asks me to set the inputs from the list of shapefiles that are in my MapCanvas. To be more specific, I always want vectorlayer_1='pst' and vectorlayer_2='astenot'.

This is the script from my model. Can anyone tell me how I can do that? I'm wondering if there is a way to do that directly in graphical modeler.

##difference_between_pst/astenot=name
##vectorlayer_1=vector
##vectorlayer_2=vector
##output_alg0=output vector
##output_alg1=output vector
##output_layer_alg2=output vector
##output_layer_alg3=output vector
outputs_0=Processing.runalg("qgis:difference", vectorlayer_1, vectorlayer_2, output_alg0)
outputs_1=Processing.runalg("qgis:difference", vectorlayer_2, vectorlayer_1, output_alg1)
outputs_2=Processing.runalg("qgis:saveselectedfeatures", outputs_0['OUTPUT'], output_layer_alg2)
outputs_3=Processing.runalg("qgis:saveselectedfeatures", outputs_1['OUTPUT'], output_layer_alg3)


texts:
I am very new to QGIS and GIS itself. I have a shapefile, which was created using Pulkovo 1995/GK zone 9(EPSG:20009) and I want to transform it to new CRS, specifically WGS84(EPSG:4326). So the main question is what is the best way to do so?

I have tried to use the Reprojection Layer feature but the coordinates of my layer doesn't change, even if my CRS changed formally.

For example: I am creating new project and define the project CRS as Pulkovo1995. Then I load the original layer and everything seems OK. In the next step I use Reprojection Layer and my layer disappears until I change the project CRS to WGS84. But the coordinates of the objects remains the same. I cannot see this layer on ESRI maps nor load it into Google Earth. Save As function doesn't help also.


texts:
I have a raster (TIFF file) and I would like to draw a mask on it.
So, I basically create a shapefile and draw polygons on the area I wanted to mask.
Then, I would like to rasterize this shapefile by creating a raster with the same extent than the raster file. This new raster must have only 0 and 1 values (1 for accepting value, 0 for not accepting).

The problem is that with the common Vector to Raster on the default plugin of QGIS, if I select the raster to get the extent (on the options), it's creating a TIFF with the initial raster value and change only the value on the polygon.

I can't draw these kind of mask on PhotoShop because the TIFF is on int16 and so, if I open it on PhotoShop I have a black image.

At this stage, the image isn't referenced. I just want to overlay a mask of 0-1 value upside the initial raster.

Is there a simple way to do tha? I can't find anything on the web... Plugin "Mask" don't do what I want.


texts: Microsoft visual studio run time library error in Qgis 1.8

texts:
When i open Qgis project, I get an error from Microsoft Visual C++ Runtime Library: Runtime Error! For other machine its working fine.Since I have uninstalled and re-installed the Qgis,but still same problem.What can I do for fixing this?

OS: Windows XP SP2

Thank you for your help.


texts: Setting layer identifiable, seachable and removable with Python in QGIS 3?

texts:
With PyQGIS, I want to set project layer properties (Project Menu > Properties > DataSource) such as Identifiable, Searchable and Required (for Read Only, the method is QgsVectorLayer.setReadOnly() (link)).

In the API doc, for the method QgsProject.nonIdentifiableLayers(), it says:

  Deprecated since version since: QGIS 3.4 use QgsMapLayer.setFlags() instead

But there is no clear documentation (for me) or example for how to use this new method QgsMapLayer.setFlags().

I would like to understand how to use this method.


texts: gIntersection results in a new dataframe with more rows

texts:
This is my first time using rgeos::gIntersection, and I have a dataset with that I rasterized using rasterFromXYZ, then turned into a spatial polygon dataframe. I wanted to intersect this SP dataframe with a shapefile (in my case different nation's administrative regions). Reading the documentation for gIntersection, the result should respect the original dimension of my rasterized spatial polygon dataframe. However for some reason, while this is true for most of the countries, the result of the generated intersection has a lot rows than the initial dataset. I could not figure out why that is and why I only see this behavior for only one country.
The rasterized spatial polygon dataframe has 11073 rows, and the shapefile has 773 features. The result of the intersection generated 19267 rows of data. I have no idea why there are now an extra ~8000 rows of data and this is currently breaking some workflow (since I need to re-bind some of the attributes original dataset onto the resultant intersection layer).
I don't have easily reproduce-able dataset/code since the data I'm working on has PII, and it only is an issue for one country. But I wanted to see if other people know why this is happening? Am I mis-understanding how the function actually works? Is there a reason why this worked for other countries but not this one?
EDIT: after some investigation on the resultant intersected SPDF, I found that certain rows of data are duplicated (but not all rows), and the duplicated rows have unique shape area and shape len fields. I'm now super confused on why the intersected spdf generated duplicated values with unique shape areas?


texts: Check get WMS/WFS capabilities for unusual server

texts:
I have an WMS(?) server that serves WMS service, I believe, but I do not understand why this server takes queer parameters that are not documented in the GeoServer WMS reference.

Could you say how to check GetCapabilities for this server?

Map link link 1,link 2


texts: If conditional for multiple raster in ArcMap10.6

texts:
I have a list of four rasters and i want to apply a lot of if then statements using Python in ArcMap 10.6. I was able get the result using raster calculator with Con statements and i tried to use the same con statements in arcpy. But i got the error message 999998 : Unexpected Error. What could be the reason for this error? Is it possible to do the same using If then statements instead of Con in arcpy and if so how? I am completely new to coding, so it would be great if someone help me fix this :)
>>> import arcpy
>>> from arcpy import env
>>> from arcpy.sa import *
>>> arcpy.CheckOutExtension("Spatial")
>>> arcpy.env.overwriteOutput = True
>>> env.workspace = "C:\Test"
>>> Raster1= Raster("Raster1")
.................
>>> outcon = Con("Raster4 " =1,3)

Thanks a lot,
/Maya


texts: Cannot get error matrix in Google Earth Engine

texts:
I'm trying to get the accuracy of my classification on sentinel-2.
I use the code below:
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7','B8','B8A','B11','B12'];

//Split sample data for training and validation
var sample = classname.randomColumn();
var split = 0.8;  // Roughly 80% training, 20% testing.
var training = sample.filter(ee.Filter.lt('random', split));
var testing = sample.filter(ee.Filter.gte('random', split));

var training = image.select(bands).sampleRegions({
  collection: training,
  properties: ['class'],
  scale: 20
});

// Create a SVM classifier with custom parameters.
var classifier = ee.Classifier.libsvm().train({
  features: training,
  classProperty: 'class',
  inputProperties: bands
});

var classified = image2.classify(classifier);

//Display classification
Map.centerObject(sample, 11);
Map.addLayer(classified,
{min: 0, max: 1, palette: ['blue','green']},
'classification');

// Get a confusion matrix representing resubstitution accuracy.
var trainAccuracy = classifier.confusionMatrix();
print('Resubstitution error matrix: ', trainAccuracy);
print('Training overall accuracy: ', trainAccuracy.accuracy());

var validation = testing;
// Classify the validation data.
var test = validation.classify(classifier,"predicted");

// Get a confusion matrix representing expected accuracy.
var testAccuracy = test.errorMatrix("actual", "predicted");
print('Validation error matrix: ', testAccuracy);
print('Validation overall accuracy: ', testAccuracy.accuracy());

When I trying to get the error matrix, GEE shows

What can I do?


texts: Selecting Feature by attribute OpenLayers

texts:
I'm attempting to limit the features on a vector layer by searching by some attribute (preferably name). I've managed to get the array of features but cannot seem to work out how to iterate over them until I find one with a specified name.

I saw in earlier versions of OpenLayers there were functions like getFeatureByAttribute which could help but there doesn't seem to be anything like this any more.


texts: Pull Data from Specific Excel Cell into ArcPy

texts:
I am new to Python and ArcPy, and wish to create a short script that pulls in a value from one specific Excel worksheet cell (namely, A2) and then performs an 'if' statement, whereby if the value of A2 is greater than 1, then a certain text string is returned, but if it is less than or equal to 1, then a different text string is returned.

So far, my code looks like this:

import arcpy
from arcpy import env
arcpy.env.workspace = r"C:\MMO_Model_Test\testdata.xlsx"
inputData = "data"

So I have defined my workspace, and set the relevant worksheet as a variable. I now need to reference the cell 'A2' as part of an 'if' statement. This is probably simple for someone who knows how, but I haven't managed to find a relevant example online. Any help is hugely appreciated.


texts: QGIS: drag zoom via keyboard shortcut/ mouse wheel?

texts:
To zoom to a rectangle in ArcGIS, I always clicked and held the mouse wheel while pushing the control key. This was a very efficiant way for me to navigate through a map. Is there anything similar in QGIS that avoids to select the tool?

EDIT:

Maybe I have to be more precise in my question, since I still have the same problem:

My initial situation is: I'm in the "select feature" mode and want to stay in this mode. Now I want to zoom in a rectangle BUT without having to click on the "zoom in" button but only with the use of a hotkey/mousewheel. In ArcGIS I always used the mousewheel-click, which allowed me to be very quick in selecting features or editing features AND navigating through the map (avoids to go to the tool bar again and again).

So: with the left mouse button, I select features (active tool in the toolbar) and with the mousewheel, clicking and holding I "panned" (not active in toolbar!) and with Ctrl + mouse-wheel-click and holding I "zoomed in rectangle".

Is there any possibility to define the hotkey in QGIS?


texts: How to convert a line shapefile to a polygon shapefile?

texts: I am using Arcgis 10.1 and have an existing line shapefile that i would like to convert to polygons.


texts: Finding documentation on ArcGIS Server's Tilemap API call?

texts:
While I was looking in Developer console, for a JavaScript App which was consuming a Tiled Map Service, I saw that it was making calls to an URL like:

https://services.arcgisonline.com/ArcGIS/rest/services/Reference/World_Boundaries_and_Places/MapServer/tilemap/4/8/8/8/8

When I tried to figure out what this is, the only place where I found some information, is in this Github Ticket: https://github.com/Esri/esri-leaflet/issues/240

According to this comment, this rest end point tells the JSAPI where there are tiles and where there are not.

That comment was back in 2014. Is any public Documentation available on this REST Endpoint?


texts: Problem creating point shapefiles programmatically using python in QGIS 2.2

texts:
I am trying to create a point shapefile in the memory buffer (for a temporary label) and then zoom to the extents on this buffer. The working (in version 2.1) test code for this is:

def originalFunction():
    vl = QgsVectorLayer("point?crs=epsg:27700", "temporary_points", "memory")
    pr = vl.dataProvider()
    # Enter editing mode
    vl.startEditing()
    count = 0
    while count < 10:
        # add a feature
        fet = QgsFeature()
        fet.setGeometry( QgsGeometry.fromPoint(QgsPoint((376000+count),(408000+count))) )
        vl.addAttribute(QgsField("mytext", QVariant.String, "", 255))
        vl.addAttribute(QgsField("mytext2", QVariant.String, "", 255))
        fet.initAttributes (2)
        fet.setAttribute(0,"hello")
        fet.setAttribute(1,"hello")
        pr.addFeatures( [ fet ] )
        count+=1
    # Commit changes
    vl.commitChanges()
    QgsMapLayerRegistry.instance().addMapLayer(vl)
    canvas = qgis.utils.iface.mapCanvas()
    qgis.utils.iface.setActiveLayer(vl)#This is done so the zoom to extent works
    qgis.utils.iface.zoomToActiveLayer()#Unless this happens the extent seems to #include the origin
    canvas.zoomByFactor(1)#this time the zoom works

This code should create a line of points diagonally and zoom to their extent in version 2.1.
When I try and run the same code in the new version (2.2 primarily but also 2.3) it does not work (no points are generated, and the fix for zooming to extent no longer seems to be working either).

I have looked at the cookbook here: http://www.qgis.org/en/docs/pyqgis_developer_cookbook/vector.html#modifying-vector-layers and tried to use that as a guide to rebuild the function. However my result (below) does not seem to be operational, despite deviating from the example code as little as I can (I found one of the functions in the example did not exist for the object it was being called on so had to replace it as shown in my comments).

New Function:

def testFunction():
    # create layer
    vl = QgsVectorLayer("point?crs=epsg:27700", "temporary_points", "memory")
    pr = vl.dataProvider()
    caps = vl.dataProvider().capabilities()

    if caps & QgsVectorDataProvider.AddAttributes:
        res = vl.dataProvider().addAttributes( [ QgsField("mytext", \
    QVariant.String), QgsField("mytext2", \
    QVariant.String) ] )
    print vl.updateFields()
    count = 0
    while count < 10:
        if caps & QgsVectorDataProvider.AddFeatures:
            feat = QgsFeature()
            feat.initAttributes (2)
#            feat.addAttribute(0,"hello")#Did not work, gave error: object has no attribute 'addAttribute'
            feat.setAttribute(0,"hello")
            feat.setAttribute(1,"hello")
            feat.setGeometry( QgsGeometry.fromPoint(QgsPoint((376000+count),(408000+count))) )
            (res, outFeats) = vl.dataProvider().addFeatures( [ feat ] )
            count+=1
    QgsMapLayerRegistry.instance().addMapLayer(vl)
    canvas = qgis.utils.iface.mapCanvas()
    qgis.utils.iface.setActiveLayer(vl)#no longer fixes the zooming to extent bug
    qgis.utils.iface.zoomToActiveLayer()
    canvas.zoomByFactor(1)
    canvas.setExtent(vl.extent())

Thanks for reading, I apologize for the length but I wanted to include code so that others can replicate the problem and see where I am going wrong!


texts: How to add additional Coordinate systems to ArcGIS Server WFS service?

texts: Has anyone added additional coordinate systems to an external WFS GetCapabilities file? I need to support webmercator to a wfs service but cannot find any documentation on how to do it. Has anyone an idea of how to do this?


texts: Why is QGIS not recognising altitude from CSV file?

texts:
I used the Google Earth method for creating a KML file of points with latitude, longitude, and altitude as necessary information. After that I converted the KML file to GPX and then to CSV, through the tcx converter. The converter returned the correct altitude expected for my region of study, but when the CSV file is opened in QGIS and used to generate a DEM through Natural Neighbor, everything goes wrong... The altitude is misread when the label is created -- the altitude should be at around 850 and is being read as 160...

How can I fix this?


texts: Replacing HTML links in attribute table Field Calculator in QGIS

texts:
I am having trouble doing a simple replace of a string in QGIS attribute table field calculator. My table has a field "tbljoined_link" that contains HTML links like Luzon. There are, however, several rows that have blank invalid/blank links like  which I want to remove in bulk.
I have tried using the following expressions but these all inexplicably replace ALL the cells in the "tbljoined_link" field (including Luzon) with NULL. I have already escaped the single quotes in the expressions below so perhaps the error has something to do with the angle brackets?
replace(  "tbljoined_link" ,'',null)

regexp_replace( "tbljoined_link" ,'',null )

regexp_replace( "tbljoined_link" ,'^$',null )


texts: Iteration of AddJoin to join CSV files and a geodatabase feature class file

texts:
I have a geodatabase feature class file and I want to join the feature class file to CSV files respectively.

For example:

Join the geodatabase feature class and CSV 1   =>  save as CSV 1 in gdb.

Join the geodatabase feature calss and CSV 2   =>  save as CSV 2 in gdb.

....
And so on until CSV 50.

Here is my code:

# Import modules and set up
import arcpy
import os

arcpy.env.overwriteOutput = True
csvfolder = "X:/SDEmanagement/imports/ACS2017/Final May 16/Test Export/"

#List csv files
csvlist=[]
for file in os.listdir(csvfolder):
    if file.endswith(".csv"):
        csvlist.append(os.path.join(csvfolder, file))

#The geodatabase feature class file to join
allgeo = "X:/SDEmanagement/imports/GEOID_joinprep.gdb/AllGeo"

#Geodatabase location to save the joined files
gdb ="X:/SDEmanagement/imports/JoinedGeoTable_topic.gdb/"

#Loop join - CSV files and All Geo file respectively.
for csv in csvlist:
    arcpy.MakeFeatureLayer_management(allgeo, os.path.splitext(os.path.basename(csv))[0]+'_lyr')
    arcpy.MakeTableView_management(csv, os.path.splitext(os.path.basename(csv))[0]+'_csv')
    arcpy.TableToGeodatabase_conversion(os.path.splitext(os.path.basename(csv))[0]+'_csv', gdb)
    arcpy.AddJoin_management(os.path.splitext(os.path.basename(csv))[0]+'_lyr','GEOID', gdb+os.path.splitext(os.path.basename(csv))[0]+'_csv', 'GEOID')
    arcpy.CopyFeatures_management(os.path.splitext(os.path.basename(csv))[0]+'_lyr', gdb+os.path.splitext(os.path.basename(csv))[0])
    arcpy.Delete_management(gdb+os.path.splitext(os.path.basename(csv))[0]+'_csv')

I got an error as below on the second or third iteration of the loop:

Error - Failed to convert Social_EdAttainment_csv. ERROR 001156: Failed on input OID -1, could not write value '13001950100' to output field GEOID
Failed to execute (CopyRows).

The error occurs in the TableToGeodatabase Function in the iteration.

Strangely, when I move the second file which gave me an error before to be in the first item in the list which makes the file to go into the loop in the first iteration, it works. So I guess the error does not come from the file itself. I just guess this is related to join in-memory files.

I tried to change personal geodatabase to file geodatabase but I got the same error.

I tried JoinField but also got the same error.


texts: Returning feature at address using ArcGIS REST Geocoding Service?

texts:
I have a geocoding service running on my ArcGIS Server Advanced 10.2.2. One of my supervisors tasked me with figuring out if there was a way to use a REST call to return a feature that is located at a certain address.

For example I have a feature service set up with a few polygons, is there a rest call that if I give it an address it would return the polygon or polygons that are at that address JSON.


texts: How to change placement for label of a polygon in mapserver with SLD

texts:
I want to change the placement of the label of special polygon by sending SLD to WMS of Mapserver. But it seems as it has no effect. Is label placement in polygons support by mapserver?
        <TextSymbolizer>
        <Font>
            <CssParameter name="font-family">arial</CssParameter>
            <CssParameter name="font-size">20</CssParameter>
        </Font>
        <LabelPlacement>
            <PointPlacement>
                <AnchorPoint>
                    <AnchorPointX>0</AnchorPointX>
                    <AnchorPointY>0</AnchorPointY>
                </AnchorPoint>
                <Displacement>
                  <DisplacementX>0</DisplacementX>
                  <DisplacementY>50</DisplacementY>
                </Displacement>
            </PointPlacement>
        </LabelPlacement>
        <Fill>
            <CssParameter name="fill">#000000</CssParameter>
        </Fill>
    </TextSymbolizer>


texts: Creating randomly shaped clumps of cells in a raster from seeds of 1 cell/pixel?

texts:
As my title says, I wish to "grow" clumps of cells from seeds within a raster. My base raster is full of 1's and 0's, 1's indicate land and 0's sea/NA areas. From the 1's I wish to select 60 random pixels/cells as my seeds, and then randomly grow a connected clump of a pre-defined no. of pixels/cells limit from that seed. I've heard the technique can be referred to as 'spread dye' but haven't had any luck finding much on it. The seed cell would be turned to a value of 2 and then the next cell chosen from surrounding 1's would be converted to 2 as well. 2's are then unavailable to be converted in future.

This thread does help a bit, as I would also be willing to do this in R as I'm familiar with reading and manipulating GIS data in R. However, what I require is a set of rules to randomly select pixels surrounding an existing clump.

If anyone has done this more basic form of cellular automata in a GIS setting then I would appreciate any advice/guidance.

Example:

I have a target of 250 cells. I randomly select one cell that has a value of 1. This is turned to a value of 2. Then, one of the neighbours of the seed cell that = 1 is turned to a 2. Then, one of the neighbours of either the cells with a 2 value is selected and turned to a 2. This would continue until a continuous shape numbering 250 cells is reached.

Edit: Further Questions

Based on whuber's great answer, I have a few questions about the code:

How do I get allocate the values of the grown cells to just a '2'
rather than variable values which represent the order in which they
were created?
I need to create 60 clumps of cells within my area of '1's.
I have devised ways to pick random starting positions but struggle
to make it all work within a loop using the expand function which you kindly
wrote. Can you suggest a way to create 60 clumps which don't clash with each other
and are contained within the same final matrix?

Edit: Further Explanation of Problem

Each clump of cells will represent a protected area of a defined size e.g. 250 cells. Each area needs to start on & grow into cells with a value of 1 as this represents land and avoid cells with a value of 0, as this represents sea. I need to iterate this 1000 times with 60 protected areas in each iteration to create a null model, showing what distributions of these areas will be by chance. For this reason, the total count of cells across all 60 areas will need to be identical in each of the 1000 iterations so they are comparable. Therefore, it's ok if the areas touch, but if there is a collision, then ideally the clump would grow in another available direction until the target of 250 is reached. Future modifications to this process will involve varying the size of each of the 60 areas and potentially the shape although I understand this would be difficult.

Once each of these 1000 protected area networks are created, they will be used as a mask against other raster data such as biodiversity measures to see (a) whether they intersect particular species ranges and (b) what % of particular species ranges these random networks of protected areas cover.

Thanks to @whuber for your help so far, I don't expect you to spend more time helping me, but thought I would try and clarify my situation as you have requested.


texts: QGIS: Fill shapefile vectors with grayscale based on elevation

texts:
I am very new to this and feel like this is probably an easy answer, but after days of searching this forum I still am unsuccessful.

I have shapefile data with elevation and I want to create a raster that is filled in (darker with more elevation)  I have tried several things (points, grids, polygons, etc) but I imagine I am missing a step in there somewhere.  The pictures show the starting shapefile, the raster lines that I have been able to create (the lines are all there, they just dont render great when zoomed at this level), and the concept I am trying to accomplish.

What am I missing?


texts:
I'm looking to download Sentinel imagery using an API (only), and as part of this process I would like to identify which tiles lay in the area I am interested in by providing the relevant bounding box coordinates. I've tried to use the OpenSearch API with URLLib in Python, as well as the sentinelsat Python module. However, in both instances I appear to have run into some sort of rate-limiting, as although I am initially able to access the relevant data, I am greeted by a 0 results page on subsequent attempts (even using the same polygon).
I have not even hit the server with that many requests (less than 10 in the space of an hour). I can't find anything about rate limiting in the Copernicus OpenSearch API docs themselves, nor the sentintelsat module documentation.
Is there rate limiting in place, if not what could be causing this issue? Is there any other way to find relevant tiles given coordinates?
I am happy to compute this myself if up-to-date data (so excluding the NTC CSV) outlining the id and footprint of each image is available.


texts: Concurrency and lock challenge in Nonversioned editing in ArcMap

texts:
I use a non version SQLServer geodatabse with geometry type for spatial data storing.To do so,I have registered all layers in SDE repository and geodatabase and ArcMap shows all layers and edit them properly.

But to apply some business rules,I have to catch ArcMap editor control in OnCreate even of Editor Object and I have to query database inside this event,so I encounter concurrency&lock problem according to this ESRI's document:

http://webhelp.esri.com/arcgisserver/9.3/java/index.htm#geodatabases/concur-1966151650.htm

even after I changed all my codes (which access geodatabase) into ArcObject model(so I use Editor::Workspace object to retrive data using ICursor and IFeatureCursor) I still receive some errors that show me database is locked while ArcMap editing session,although ESRI says isolation level is UNCOMMITTED READ.If isolation level is uncommitted read,I can at least read data from geodatabase during the time after ArcMap creates a feature and before user commits transaction.But event a simple SELECT query is locked by ArcMap.I wonder if I can't understand right meaning of UNCOMMITTED READ isolation level or there is a problem some where else.

Any suggestion or reference to deal with this problem?

I use ArcGIS Desktop 10.2 for working with data and use ArcSDE 9.3 to create and register geodatabase.


texts:
I am trying to get a where_clause working on a UpdateCursor (arcpy.da). I can get it to work fine with a SearchCursor, but for some reason the same clause throws a RuntimeError: General function failure error.

import arcpy
from arcpy import env
env.workspace = "C:/mydata/mylocation/Results"
fc = "airports.shp"
delimfield = arcpy.AddFieldDelimiters(fc, "STATE")
queryString = delimfield + " <> 'AK'"
cursor = arcpy.da.UpdateCursor(fc,["STATE"], queryString)
for row in cursor:
    print row
    row[0]="AK"
    cursor.updateRow(row)
del row
del cursor


texts: Getting table with coordinates and pixel value for polygon from Google Earth Engine

texts:
I would like to have a table or a matrix with the coordinates and value of each pixel on a given polygon. I am working with the Global Forest Change data specifically with the 'loss' layer. In this layer, each pixel has a value of 1 if the forest has been cleared or 0 if it was not.
What I have done so far makes no sense as I am a beginner. The code below is the one offered by the Google Earth Engine Tutorial.
var gfc2018 = ee.Image("UMD/hansen/global_forest_change_2018_v1_6");

Map.addLayer(gfc2018.mask(gfc2018), {
  bands: ['treecover2000'],
  palette: ['000000', '00FF00'],
  max: 100
}, 'forest cover masked');

var treeCover = gfc2018.select(['treecover2000']);
var lossImage = gfc2018.select(['loss']);
var gainImage = gfc2018.select(['gain']);

// Use the and() method to create the lossAndGain image.
var gainAndLoss = gainImage.and(lossImage);

// Add the tree cover layer in green.
Map.addLayer(treeCover.updateMask(treeCover),
    {palette: ['000000', '00FF00'], max: 100}, 'Forest Cover');

// Add the loss layer in red.
Map.addLayer(lossImage.updateMask(lossImage),
            {palette: ['FF0000']}, 'Loss');

// Add the gain layer in blue.
Map.addLayer(gainImage.updateMask(gainImage),
            {palette: ['0000FF']}, 'Gain');

// Show the loss and gain image.
Map.addLayer(gainAndLoss.updateMask(gainAndLoss),
    {palette: 'FF00FF'}, 'Gain and Loss');

// Create a polygon
var region = ee.Geometry.Rectangle(-75.2806, 3.1209, -75.0554, 3.2413);

// Extract pixel data from polygon
...?```


texts:
I've downloaded monthly SST data (11u, nighttime) for July 2014.

https://oceancolor.gsfc.nasa.gov/cgi/l3?per=MO&prd=NSST_sst.nc&sen=A&res=4km&num=24&ctg=Standard&date=1Jul2014

The data is stored as a NetCDF (.nc) file. So, I've used the following code to import the file into R.

library(ncdf4)
library(raster)
sst<-raster(path.expand("~/Remote Sensing/SST/SST_July2014_L3_MO_4km.nc"))

Next, I checked the file.

sst

This was the output:

class       : RasterLayer
dimensions  : 4320, 8640, 37324800  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -90.00001, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : ~\Remote Sensing\SST\SST_25June2016_L3_8D_4km.nc
names       : Sea.Surface.Temperature
zvar        : sst

Does the raster() function in ncdf4 automatically project the file? I'm confused because the MODIS SST metadata states that the projection is Equidistant Cylindrical.
Sorry if this question is very elementary, but this is my first time working with rasters or MODIS data.

EDIT

I successfully specified the projection using the code provided by @Phil

sst<-raster(path.expand("~/Remote Sensing/SST/SST_July2014_L3_MO_4km.nc"))
sst@crs <- sp::CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs")

My other data files are in latlong WGS84, so I reprojected the SST raster.

sst<-projectRaster(sst, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

I then tried to crop the bathymetry raster by the extent of my polygon layer called 'polygonsp.' (code for polygonsp below)

library(sp)
vertices<-c(-131.8246,-131.8246,-127.6479,-127.6479,50.14335,54.53591,54.53591,50.14335)
vertices<-matrix(vertices,4,2)
polygon<-Polygon(vertices)
polygono<-Polygons(list(polygon),1)
polygonsp<-SpatialPolygons(list(polygono))
proj4string(polygonsp)<-CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

sstcrop<-crop(sst,polygonsp)

However, I now get an error that says "Error in .local(x, y, ...) : extents do not overlap".

This code worked for me before I added the line specifying the projection of the SST data as equidistant cylindrical. Any idea what is happening?


texts: Rule-based rendering in qgis2web client is not working

texts:
I made a simple map with polygons of the staates of africa and used the rule-based renderer for different colors to tint them. When I tried to load or export this map with qqis2web I didn't get displayed any colors but only the contours. I tried several very simple rules but in every case the result was always the same. If I used the categorized renderer for that everything worked fine.
Isn't the qqis2web able to display rule- based rendering or what could be the matter in that case?


texts: Vector points from CSV file wrong location in QGIS

texts:
I am using QGIS 1.8 (I know old version), and uploaded a CSV file with coordinates. The underlying base map is from Google, using the OpenLayers plugin. This base map uses WGS84-EPSG:3857, and my vector points from the CSV are using WGS84-EPSG:4326. This difference shouldn't be a problem I thought to have read in here:

EPSG 3857 or 4326 for GoogleMaps, OpenStreetMap and Leaflet

The strange thing is that, my points appear on the map just about 30 kilometres next to the actual location. If I used the wrong projection or coordinate system then the point would have been on the other side of the world. A colleague of me uses the exact same file and we think he enters the same settings, at least the same CRS, and his points are in the right spot. Does anyone knows what the problem might be?


texts: PostGIS not converting GIS data into latitude/longitude values

texts:
I am running PostgreSQL and PostGIS. I am trying to do something mind-numbingly simple: determine the latitude and the longitude of a point.

Unfortunately, PostGIS stubbornly refuses to transform the SRID of the GIS data from SRID=4269 to SRID=4326, which is necessary to get the latitude and longitude of the point in question:

select pid, landdistri, ST_AsText(ST_Centroid(geom)) from "parcels_new" where pid='016-626-331';
     pid     |       landdistri        |                st_astext
-------------+-------------------------+------------------------------------------
 014-626-331 | Similkameen Div of Yale | POINT(1461547.22664562 504522.034305431)
(1 row)

select pid, landdistri, ST_SRID(geom) from "parcels_new" where pid='016-626-331';
     pid     |       landdistri        | st_srid
-------------+-------------------------+---------
 014-626-331 | Similkameen Div of Yale |    4269
(1 row)

select pid, landdistri, ST_SRID(ST_Transform(geom, 4326)) from "parcels_new" where pid='016-626-331';
     pid     |       landdistri        | st_srid
-------------+-------------------------+---------
 014-626-331 | Similkameen Div of Yale |    4326
(1 row)

select pid, landdistri, ST_AsText(ST_Centroid(ST_Transform(geom, 4326))) from "parcels_new" where pid='016-626-331';
     pid     |       landdistri        |                st_astext
-------------+-------------------------+------------------------------------------
 014-626-331 | Similkameen Div of Yale | POINT(1461547.22664562 504522.034305431)
(1 row)

select pid, ST_AsEWKT(ST_Centroid(ST_Transform(geom, 4326))) from "parcels_new" where pid='016-626-331';
     pid     |                     st_asewkt
-------------+----------------------------------------------------
 016-626-331 | SRID=4326;POINT(1461547.22664562 504522.034305431)
(1 row)

select f_table_name as name, srid from geometry_columns;
    name     | srid
-------------+------
 parcels_new | 4269
 Pins        |    0
(2 rows)

I don't get any error messages.
I am at a loss to understand why ST_Transform REFUSES TO TRANSFORM my data....


texts: polygon popupContent did not get updated

texts:
I have a $.getJSON function, which gives me a list of polygons and multipolygons, I am binding popups to each polygon like this

[inside the polygon loop]:

           var onEachFeature = function(feature, layer) {
              if (feature.properties && feature.properties.popupContent) {
                 layer.bindPopup(feature.properties.popupContent);
              }
           }

            var polyStyle = {
                    "color": "#f2f2f2",,
                    "weight": 2,
                    "opacity": 0.5,
                    "fillColor": "#f2f2f2",
                    "fillOpacity": 0.4,
                    "radius": 8,
                    onEachFeature: onEachFeature
            };
            var geojsonFeature = {
                    "type": "Feature",
                    "properties": {
                        "name": feature_name,
                        "popupContent": feature_content //displays uniq data belongs to the polygon
                    },
                    "geometry": polyData
            };

    geojsonLayer = L.geoJson(geojsonFeature, polyStyle).addTo(map);

The popup contents on all the polygons display the same content, which is the content of the last polygon data in the loop. This is wrong, each popup content should show different data.


texts:
I want to add layout number to each layout in the print composer.  I've cant find a variable with that information in it that I can access from the expression builder.

Generally our projects are built with many single page layouts.  They need a unique ID per project, but this ID can overlap with other projects.  So I figured we could automate away the task of making sure each layout has a separate ID number printed on it, only I can't figure a way to get at the data within the expression builder.

Does anyone have an idea on a way to get a unique ID onto each print layout?  (would prefer an incriminating one like 0,1,2,3,4)


texts: Spatial Join Points and Polygons in ArcGIS Desktop?

texts:
I am trying to create a grid road map with an index table for street look up.

I am using this tutorial by ESRI:
http://blogs.esri.com/esri/arcgis/2007/11/08/creating-street-name-indexes/

I am at the step whereby I am supposed to spatially join a point file (centroid point for each label) to a grid (which has the grid ID):

  The next step is to use the Spatial Join tool (one to many is the join
  option) and join the annotation to each grid cell they fall within.
  However, you will likely not want to work directly with the annotation
  because the spatial join options for “contains” or “within” will not
  include any annotation that overlaps more than one grid cell. Further,
  the “intersects” option will create duplicate entries in the same
  circumstances. To get around that convert the annotation to point
  features using the Feature to Point tool (Important note: I had to
  delete the Element field; the type is blob, the field isn’t needed,
  but if it’s there, the tool won’t run).

It seems to be telling me to spatially join the points to the grid in order to match up each label point with its grid ID value. However, whenever I spatially join, the output table is empty.

Any advice? I simply need to find a way to extract the grid ID of where the points falls within the grid.

I've included a photo of the spatial join parameters. After the process runs, it simply spits out a new shapefile with no features and an empty attribute table.


texts: How to Create a Custom-shaped Fishnet Grid in QGIS

texts:
I am looking to automate the create of a fishnet grid over a large area (about the size of Rhode Island) for a mapping study. I want to use hexagonal grid cells, but the kicker is I need each of these hexagonal grids again broken down into their 6 hexants. The center convergent vertex will be a setup location with mapping done out to the extent of each hexant. Sort of like flying the petals from the pistil of a flower.

I don't see custom grid shapes in QGIS. I'm wondering how I can create this grid layout.
My first thought was to lay down hex grids and then do a series of layer copying with systematic offsets and a final merge to get what I need. Or to convert grids to polylines, copy those lines with the bearings I want, and repaste/offset those. But there's got to be other tools or a simpler process out there.


texts:
I am experimenting with Potree and I am very pleased with the viewer. However, I would like to customise the sidebar but seem to be having trouble. The source code in the developer tools seems to reference files in an src folder which don't seem to be referenced anywhere I can get at them to edit. I do have access to the sidebar.html but I'm not sure what this does or how it's referenced if at all.
Is there an obvious way I can use the code from the sidebar.html in the webpage generated when the pointcloud was converted?


texts:
I've got many GPX files in the following format:

<?xml version="1.0" encoding="utf-8"?>
<gpx xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xsd="http://www.w3.org/2001/XMLSchema" version="1.1" creator="c174a4ce-d611-43ef-b4d9-e4122bc8f2a4" xmlns="http://www.topografix.com/GPX/1/1">
  <rte>
    <name>Withens Clough, Cragg Vale and Stoodley Pike</name>
    <rtept lat="53.705918666306154" lon="-2.02137437495311">
      <ele>268</ele>
    </rtept>
    <rtept lat="53.705457119907116" lon="-2.0253609418514311">
      <ele>293.8</ele>
    </rtept>
    ....
  </rte>
</gpx>

I use this is a website to generate a route using OS Openspace maps. However, this doesn't create an elevation profile.

Is there a web service or something I can run on my web server that can generate GPX elevation graphs for embedding into webpages? My web server runs PHP.

I've found plenty of Google map ones but don't need the map, just the graph, but in the implementations I've looked at it seems like one depends on the other and don't want to compromise loading times by including two maps or remove the OS maps which are a lot better, in my view.


texts: Copy features between file geodatabase and ArcSde by using ArcObjects

texts:
I want to transfer all contents of a table from file geodatabase (*.gdb) to a table in ArcSde. In fact the data is meant to be replaced, so I need also to clear the table in ArcSde. Both tables have same structure. The target table is used in a map service.

I started with ConvertFeatureClass and the effects were pretty good - it was working, but each time it was creating a new table. Every try to append the data to already existing table ended with errors. Here's my code - http://pastebin.com/RWhMTe4p

Since this was not doing what I wanted, I tried something different. This time, I used IGeoDBDataTransfer interface. My efforts ended up with an error being thrown on the final Transfer method. That's my code - http://pastebin.com/MG3KUHdS The COMException is: "Error -2147216061: Exception from HRESULT: 0x80041543

Can someone please explain how to do this task?


texts: Tiled map service base layer does not display properly in ESRI Javascript map

texts:
I am adding a map to a website built in asp.net and Javascript. When I add the esri.Map using a tiled map service layer base map, the map overlaps the map div and what I assume are tiles show up out of order and with spaces between them. Picture a large checkerboard with transparent squares alternating with map squares and Antarctica in the map square to the right of North America and Greenland at the bottom. There also appears to be some metadata text appearing in the div in a transparent square.

There is another overlapping div that is turned off and the map div is turned on in response to a button click to display the map.

A dynamic map service layer displays properly, but it behaves a bit oddly on zoom, extending to take up the whole page during the zoom, then snapping to the div when the zoom is complete.

The map div is being set in css to relative position, and I've tried also giving it a set size in pixels, with no change in behavior. Originally the div was not given any position or size attributes, but the behavior was the same. I get the same result in Opera, Chrome, and Firefox. Nothing shows up in Explorer, but I think that may be a security setting issue.

I'll  try to attach code and screenshot of the map.  Below is the javascript code adding the map. I'll just note that I've tried various extents and spatial references, including copying and pasting from standard "add a map" tutorials, so this is just what's there now. They all look the same. I've also added the basemap with addLayer by url with the same results.
Javascript function to add map:

function MapIt(msg) {

//display map box
$("div.selection-criteria").hide();
$("div.map").show();

map = new esri.Map("map", {
    extent: new esri.geometry.Extent({ xmin: 38.705772, ymin: -84.820221, xmax: 39.610317, ymax: -83.683907, spatialReference:{wkid:4326} }),
    zoom: 3,
    basemap: "streets"
});

}

ASP code for the content and divs, with the fieldsets for the div that gets hidden when the map is show obfuscated:

<asp:Content ID="Content1" ContentPlaceHolderID="contentMain" runat="server">

<fieldset class="somefields">
</fieldset>

<fieldset class="someotherfields">
</fieldset>

<input type="button" id="btnMap" class="map-button" value="Map it!" />

CSS div styles:

div#map{position:relative; width:800px; height:400px;}

![map]: http://i.imgur.com/65yxznk.jpg


texts: Exporting SQL Server .bak files to a shapefile?

texts:
I have no experience using Microsoft SQL Server at all. Someone gave me a .bak file of their database and said that I can extract their spatial data that way. They said that they would not export them as JSON, KML, or shapefiles for me.

So... how can I open the .bak file and export the tables/data that I need as shapefiles? I am currently downloading SQL Server Express. I have very limited knowledge of SQL.


texts: Area calculation of an object without coordinates

texts: I have a photo of an object with a scale and I want to calculate its area. I guess I have to create a polygon of the area to be calculated to use the 'field calculator'. But how do I tell the program what a centimeter (cm) is?


texts: Calculating imprecision introduced by using WGS 84 coordinates as x,y in raster?

texts:
I have several coordinates provided in WGS 84 which I need to rasterize. For that purpose I simply interpret the lat,lon coordinates as x,y coordinates and I apply some scaling. So e.g. anything in the rectangle described by (54.3460005,10.1494873)(54.2984447,10.2002991) will be scaled to the rectangle (0,0)(10000,10000). In further processing steps this raster is translated back into the WGS 84 coordinate system.

I am aware of that this simple translation would be problematic for larger areas but I only look at smaller urban areas so I think the error introduced is acceptable. Further, visual inspection showed that the difference between my raster method translated back into WGS 84 coordinates and the original coordinates is practically neglectable. Yet I feel uncomfortable doing so without knowing the exact error I introduce by this method. At least I want to report it.

How can I calculate on a theoretic basis how large the error is which I introduce by simply interpreting lat,lon (WGS 84) values as x,y values?

EDIT: The error I am concerned about is the level of distortion of a pixel. While in the northern part of the map a pixel covers a smaller area (when being in the Northern hemisphere like my example), in the southern part the area covered is larger.


texts: Multi-page export in *.tiff and other image formats using Export Map?

texts:
Is there a way to export multiple maps in *.tiff using "Export Map" and setting Data Driven Pages?

Because Data Drive Page work in "Export map" only with PDF.


texts: ogr2ogr error importing shapefile into PostGIS : "numeric field overflow"

texts:
I'm trying to import a shapefile into a PostGIS database using ogr2ogr as follows:

ogr2ogr -lco GEOMETRY_NAME=geom -lco LAUNDER=NO -a_srs "EPSG:4326" -f "PostgreSQL" PG:"dbname=db" test.shp

This is after trying the import with shp2pgsql and encountering an issue OGR integer type to PostgreSQL BIGINT?.
The shapefile was a Tiger/Line file circa 2000 if I remember correctly. But I am getting the following error:

ERROR 1: COPY statement failed.
ERROR:  numeric field overflow
DETAIL:  A field with precision 19, scale 11 must round to an absolute value less than 10^8.
CONTEXT:  COPY test, line 1, column TLID: "142691303.00000000000"

The TLID field in question shows as double real 18 11 for type, type name, length and precision in QGIS 3/2.99. I did some quick search, and TLID seems to be TIGER/Line ID. I don't why it is shown as a double since it's an ID, or why is it causing this error for ogr2ogr.

This is tested under Ubuntu 16.04, and ogr2ogr --version shows:

GDAL 2.2.1, released 2017/06/23

Does anyone know what went wrong and how to fix this to allow importing into PostGIS?

-- EDIT --

The dbfdump -m output for the file is:

Record: 0
STATEFP: 06
COUNTYFP: 037
COUNTYNS: (NULL)
TLID: 142691303.00000000000
TFIDL: 219999921.00000000000
TFIDR: 219999382.00000000000
MTFCC: S1400
FULLNAME: Fashion Ave
SMID: 2361
LFROMADD: (NULL)
LTOADD: (NULL)
RFROMADD: (NULL)
RTOADD: (NULL)
ZIPL: (NULL)
ZIPR: (NULL)
FEATCAT: (NULL)
HYDROFLG: N
RAILFLG: N
ROADFLG: Y
OLFFLG: N
PASSFLG: (NULL)
DIVROAD: N
EXTTYP: N
TTYP: (NULL)
DECKEDROAD: N
ARTPATH: (NULL)
Shape_Leng:      42.48449176210

Record: 1
STATEFP: 06
COUNTYFP: 037
COUNTYNS: (NULL)
TLID: 141692965.00000000000
TFIDL: 219999899.00000000000
TFIDR: 219999383.00000000000
MTFCC: S1630
FULLNAME: (NULL)
SMID: 2361


texts: Can't load layers into ArcGIS Viewer for Silverlight

texts:
I am trying to load layers into ArcGIS Viewer for Silverlight. I have uploaded my files into the server successfully. However when I browse for the files I get connection error. I know the files are uploaded into the server because I can see them when I use the same file path to the server for Flex Viewer.

Thanks for the help.


texts: Geoserver: WMS layer is shifted at certain zoom scales in openlayers

texts:
We use geoserver to serve WMS tiles inside a openlayers environment. The weird thing I see is that, in our production environment, depending on the scale, data is shifted (translated) against the OSM/Google/Bing background.

See following screenshot, here you can see the pipelines positioned correctly on the plant site:

Now if I zoom out, I get the following view:

obviously wrong. I am a bit puzzled as to what could cause this.

I develop on a mac, using geoserver 2.6.3 and (unfortunately) Oracle database backend. On my mac I do no have this problem at all. So it is imho not related to my OpenLayers/javascript code.

Our production environment is a Windows 2012R2 server, with Oracle 11g2. I first ran Geoserver 2.6.3, just upgraded to 2.7.1, hoping that would remove the problem (but it did not --a lot faster though, which is nice).

I read on the geoserver site that Windows+oracle is an obscure platform, but I am kind of hoping somebody could think of something to look at that could cause this?

Things I tried/verified:

the data is the same on both platforms
the data is stored in SRS-ID 3035, and thus transformed to 900913
the cache is empty (originally the srs-id was incorrectly deduced from the table)
used the same browser
I checked the featuretype.xml on both platforms and it looked identical

At the moment I am bit stuck, my guess it is the combination of oracle and windows that causes this error? We use postgis+geoserver 2.3 on windows for another client without issues.

But I have no idea where to look for a probable cause, and eventually, hopefully a fix.


texts: Seeking Icon library for custom GIS software?

texts:
I am building a desktop application and currently utilizing a mix of icons for pan, select, zoom, attribute table, etc that have been randomly taken from the internet and modified.

The overall look is very messy and not that intuitive to our non-technical users so
I am hoping to find a set of icons that have a homogeneous look to them.

My team has a small budget to spend on buying a set of icons but, if they exist, we would prefer to use free open source icons.

Can you point us to resources for either paid or free icon options?


texts: GeoPandas imports null attribute values as zero from shapefile

texts:
I would read a shapefile as GeoPandas DataFrame using:
import geopandas as gdp
df = gdp.read_file('file.shp')

When I read the file, the null values have been loaded as zero. I need to see and distinguish null values and the real zeros.


texts:
Is there an OSM label layer available? This same question was asked 7 years ago. I neither see such a layer in the OpenLayers plugin nor QuickMapServices.
I'm working on a project with a 50% transparent layer over OSM monochrome. It would look better to have the OSM labels over my layer, with the basemap free of labels.


texts: adding barriers is not working in arcgis sdk for android

texts:
I have a program that displays this path when I give it two points shown by the red marker.

I wanted to create a barrier on that road so I added this code

        Polygon po=new Polygon();
        Point p1=new Point(4102363.9071817156,857772.619643);
        Point p2=new Point(4102416.5557467686,857799.2401230272);
        Point p3=new Point(4102436.1589281,857744.31787604);
        po.startPath(p1);
        po.lineTo(p2);
        po.lineTo(p3);
        Graphic ba_gr=new Graphic(po,null);
        NAFeaturesAsFeature na=new NAFeaturesAsFeature();
        na.addFeature(ba_gr);
        params.setPolygonBarriers(na);

params is my RouteParameters

On my previous code but my program still shows be the previous path. Did I miss something when creating the barriers?

The polygon is in this position,

I tried it on another program.


texts: caching part of a postgis layer in mapserver

texts: I'm creating a layer in mapserver using postgis. The layer consists of many polygons. Sometimes some of those polygons will change. Is it possible to cache the layer tiles using something like tilecache and only regenerate tiles containing the changed polygons? Or do I need to regenerate the whole layer?


texts:
I have nc file where time unit is expressed as year and month in format "month as %Y%m.%f", type: character. For example, 199604.5 represents the year 1996, month 04. How to convert time in date format?

I tried next:

as.Date(ncin$dim$time$vals, "%Y%m.%f")

But I got the next error:

Error in charToDate(x) :
  character string is not in a standard unambiguous format

I also tried next:

as.Date(ncin$dim$time$vals)

but I got the wrong dates.

"2516-07-01" "2516-07-02" "2516-07-03" "2516-07-04" "2516-07-05" "2516-07-06"


texts:
I am trying to understand the process of Radiometric Calibration of multi-spectral images acquired using Parrot Sequoia sensor with firmware v1.2.1.

Given a radiometric calibration target and the Albedo values for all bands of interest (Red, NIR, Red-Edge), how to apply this information to get a radio-metrically calibrated image.
What is the formula to apply radiometric calibration?


texts: Is it permissible to update features after the feature cursor has been released?

texts:
In ArcObjects, is it permissible to modify and update a feature after the feature cursor (which was used to retrieve the feature) has been released?

That is, do features in any way retain a link to the cursor through which they were retrieved?

// retrieve a feature through a recycling feature cursor:
IFeatureClass featureClass = …;
IFeatureCursor featureCursor = featureClass.Search(nothing, true);
IFeature feature = featureCursor.NextFeature();

// release the feature cursor:
Marshal.FinalReleaseComObject(featureCursor);

// QUESTION: can the feature still be modified and updated?
feature.Value(…) = …;
feature.Store();

I'm programming against ArcGIS versions 9.3 and 10.
Does it possibly make a difference whether the cursor is recycling or non-recycling?

P.S.: I have tried this under both ArcGIS versions 9.3 and 10, and it seems to work fine; however, that doesn't mean that this will always work, nor that it should be done. I haven't found any hints about this issue in ESRI's online API documentation.


texts:
When I issue "CREATE EXTENSION postgis", I'm getting :

  curs.execute('CREATE EXTENSION postgis')
  psycopg2.InternalError: could not load library "/usr/lib/postgresql/9.5/lib/postgis-2.3.so": /usr/lib/liblwgeom-2.3.so.0: undefined symbol: GEOSClipByRect**

I use psql (PostgreSQL) 9.6.1. Does anybody know how can I fix it?

I use following comments for installing postgis.

sudo sh -c 'echo "deb http://apt.postgresql.org/pub/repos/apt xenial-pgdg main" >> /etc/apt/sources.list'

wget --quiet -O - http://apt.postgresql.org/pub/repos/apt/ACCC4CF8.asc | sudo apt-key add -

sudo apt-get update

sudo apt-get install postgresql-9.5-postgis-2.2 pgadmin3 postgresql-contrib-9.5


texts: How to convert the DN values of Resourcesat-2 LISS-3 images to Radiance?

texts:
Does somebody knows How to convert to Radiance the Resourcesat-2 LISS-3 images?

I found a paper "Pre-flight Radiometric and Spectral Calibration of Resourcesat-2ALISS3* payload", where they determined some parameters. But I´m not sure. I´ve also tryed the SACRS2, of Vishal Pathak called at the paper "Retrieval of Surface reflectance from Resourcesat-2 AWiFS, LISS-3 and LISS-4 data using SACRS2 scheme", but I was not able to put the right input data (RPC, META, geom?)


texts: Get spatialReference Domain, with XY Precision

texts:
I'm trying to get a spatialReference domain from prj file by using the function :

private void ConstructCoordinateSystem(bool highPrecision)
{
    // SpatialReferenceEnvironment is a singleton object and needs to use the Activator class.
    Type t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
    System.Object obj = Activator.CreateInstance(t);
    ISpatialReferenceFactory3 spatialReferenceFactory = obj as
        ISpatialReferenceFactory3;

    // Create a geographic coordinate system from a .prj file that contains the well-known string.
    // You might need to update the file location.
    ISpatialReference3 spatialReference =
        spatialReferenceFactory.CreateESRISpatialReferenceFromPRJFile(
        "C:\\Program Files\\ArcGIS\\Desktop10.0\\Coordinate Systems\\Geographic Coordinate Systems\\World\\WGS 1984.prj")as ISpatialReference3;

    // Determines whether you are constructing a high- or low-precision spatial reference.
    IControlPrecision2 controlPrecision = spatialReference as IControlPrecision2;
    controlPrecision.IsHighPrecision = highPrecision;

    ISpatialReferenceResolution spatialReferenceResolution = spatialReference as
        ISpatialReferenceResolution;
    ISpatialReferenceTolerance spatialReferenceTolerance = spatialReference as
        ISpatialReferenceTolerance;

    // The following three methods are key: construct horizon, then set the default x,y resolution and tolerance.
    spatialReferenceResolution.ConstructFromHorizon();
    spatialReferenceResolution.SetDefaultXYResolution();
    spatialReferenceTolerance.SetDefaultXYTolerance();

    // Check the x,y spatial domain.
    double xMin;
    double xMax;
    double yMin;
    double yMax;
    spatialReference.GetDomain(out xMin, out xMax, out yMin, out yMax);

    //System.Windows.Forms.MessageBox.Show("Domain : "+ xMin +", "+  xMax + ", "+ yMin +", "+ yMax);
}

When I use a prj file other than WGS, it works perfectly.

But when I use a WGS file, in this case I need to add a test hasXYPrecision before calling spatialReference.GetDomain().

The problem is that the property hasXYPrecision returns always false, in both cases (WGS file and non-WGS file)

Any ideas?


texts:
I have tables/layers with hierarchical data in it with and without geo-positions.
Example without:
fid, projectnames <table1>
 -> fk_table1_fid, projectpart, display ... <table2>

Example with:
fid, cable, fk_table2 <table3>

I know I can edit the table1,2,3 by double clicking on a entry in the tableview or by clicking on the feature in map-view.
But for tab 1+2 I do want to have a better workflow/handling. I dream of a GUI where I can add a project and then add parts to it without closing the GUI, open table double click or add... The most important thing I want to have is a much simpler way to use the display-field.
With this field I control the displaying of features in my map. I do want to change it individually on parts-table and for all parts from table1's projects. (If I click to hide a full project, then all parts have to become false too).
Is this reachable without writing a plugin? I am a noob in Python... I know a little bit about C++ but more with PHP or VB.net.
I tried to make a plugin with plugin builder but it can not found pyrcc5, I installed it resp. There are some .bat with this names on my system, one of this is in %PATH% but it can not be found... at this point I am lost...
Edit: I do want to start editing from a button in a toolbar, or mostly welcome from a dockable GUI/Tab.


texts: Road Shield Labels not displaying properly

texts:
I am operating in ArcGIS Desktop 10.3. I have created three different road shield labels for the three main road types; interstate highway, us highway, state highway. However, the orientation of the symbols is not desirable.

I want the road shields to be straight like this:

The strange this is, although these are two different map projects (MXDs), the shapefile is exactly the same and as far as I can tell all the label placement settings are set the same way.

Any ideas?


texts:
Just wondering if there is a method of creating a layer for each classification range in a TIN.

Another way to state my problem:
I would like to make each colour code a layer that I can "turn on and off" i.e. it will have a check box at the front.

Thanks
David


texts: Plotting shapefile on map using ArcGIS API for JavaScript?

texts: We are using Dojo framework to develop our portal. We have bound the ArcGIS map in our web portal. Is there any ArcGIS JavaScript API which will accept the path of the shapefile and plot that shapefile directly on map. If it is free API then that will be great and if it is paid API then how can we purchase that?


texts: Converted shapefile to GeoJSON - how to put this on OSM map?

texts:
I converted a shapefile to a GeoJSON file but I am having trouble adding this data to a map using Leaflet.

The osmscarb.js:

var osmscarb = {
"type": "FeatureCollection",

"features": [
{ "type": "Feature", "properties": { "Id": 0, "POINT_X": -8826927.77461, "POINT_Y": 5417539.05834, "INTXN": "Luttrell\/Danforth", "FILENAME": "Luttrell Loop.jpg" }, "geometry": { "type": "Point", "coordinates": [ -8826927.774612015113235, 5417539.058339669369161 ] } }

I'm putting this in my map:

var myLayer = L.geoJson().addTo(map);
myLayer.addData(osmscarb);


texts: Create fixed number of Points within polygon

texts:
I have a polygon shapefile containing 50 features. I have another list where those 50 features are listed with how many number of points need to be created within those polygon. Suppose I have a shapefile containing A, B, C, D polygon and an excel file like following

+-----+-----+
|  A  |  2  |
+-----+-----+
|  B  |  6  |
+-----+-----+
|  C  |  9  |
+-----+-----+
|  D  |  3  |
+-----+-----+

Is there any easiest way to create those number of points within assigned polygon?


texts: how to convert back from pixel coordinates to x,y coordinates

texts:
Here is what I have done to project shapely points to pixel coordinates :

#shapely list of points
points= ['POINT (656822.0796426814 -185003.7243437266)', 'POINT
(656917.7545657885 -184985.6672704536)', 'POINT (656997.7888963958
-185001.578684116)']

#converting shapely list to np array
coords_list = [[geom.xy[0][0], geom.xy[1][0]] for geom in points ]
coords = np.array(coords_list)

# Reproject shapefile points to numpy coordinate system
proj_im = gdal.Open(input.tif')
geoinformation = proj_im.GetGeoTransform()

xorigin = geoinformation[0]
yorigin = geoinformation[3]
pix_size = abs(geoinformation[1])

coords_unp = np.roll(np.column_stack(((coords[:, 0] - xorigin) / pix_size, \
                                      (yorigin - coords[:, 1])/pix_size)),1, 1)
print(coord_unp)
>> [[ 2359.02892404  1537.68823732]
   [ 2085.43690476  2987.30828439]
   [ 2326.51892994  4199.94965723]]

How do I convert back from pixel coordinates to lat,lon?


texts: What is the best way to protect data online?

texts: What techniques do people use when protecting data (e.g., preventing a user downloading the whole dataset) they are overlaying on top of Google/Yahoo/Bing etc.?


texts:
It is possible to add a new feature to qgis web client?
I would like to add jump/zoom to selected layer. I tried to add something like this using ZoomToExtent feature with fixed coordinates, but this didn't work. I think I could get extent coordinates of selected layer but I can't pass it to ZoomToExtent. Guidelines would be nice.


texts:
I'm trying to develop a robust way to generate a UI incorporating map layer controls such as a base layer selector. I've read several posts including this recent one but I'm still not able to change the visibility of the layers and always seem to return an error 'layer.setVisible is not a function'. Once I have a reliable method, I will be developing other controls for layer-styling etc.

I've followed several tutorials and generated a layer tree with opacity controls for example, but even though all of my layers are in groups, the layer tree generates several 'undefined' layers which I can't locate. Hence the reason I would like to build this from scratch, so that I can understand how to include or exclude layers and ensure that the ui is bound to the map object.

I'm aware of Matt Walker's LayerSwitcher but again, it doesn't really do what I want and I'm struggling to understand the 'layer walking' process adequately. Below is a sample of my code:

function getBaseGroup() {
var layers = map.getLayers();
var len = layers.getLength(), l;
for (var i = 0; i < len; i++) {
    l = layers.item(i);
    var lt = l.get('title');
    // check for layers within groups
    if (lt == 'Base Layers') { // Title of Group
        if (l.getLayers) {
            var innerLayers = l.getLayers().getArray();
            var base = [];
            for (var j = 0; j < innerLayers.length; j++) {
                base.push(innerLayers[j]['B']);
            }
            return base;
        }
    }
}
}

function showBase() {
var this_ = this;
var sel = $('#baseMapLayer');
var layers = getBaseGroup();

layers.forEach(function (layer) {
    sel.append($("<option>").attr('value', layer.title).text(layer.title));
});

$(sel).change(function () {
    var selectedValue = $(this).find(":selected").val();
    console.log("the value you selected: " + selectedValue);
    layers.forEach(function (layer) {
        console.log(layer); // returns layer object
        if (layer.title == selectedValue) {
            console.log(layer.title); // returns title
            layer.setVisible(true);
        }
        else {
            layer.setVisible(false);
        }
    });
});


texts: Creating proper GeoJSON from PostGIS data using node.js

texts:
I am trying to build a GeoJSON object from a SQL query to some GIS point data in a postgis postgresql database. A snippet of my node.js app.js is below.

As it stands I understand building the type and features, but don't know how to attach a properties properties to each GeoJSON record (in the below, it all renders at the end, separate from (not collated with) the features).

THE QUESTION: What do I need to do so that the properties attach (collate) for each record in the loop that builds the GeoJSON So that it looks more like this http://www.geojson.org/geojson-spec.html#examples?

function GrabData(bounds, res){

  pg.connect(conn, function(err, client){

  var moisql = 'SELECT ttl, (ST_AsGeoJSON(the_geom)) as locale from cpag;'

  client.query(moisql, function(err, result){
    var featureCollection = new FeatureCollection();

    for(i=0; i<result.rows.length; i++){
      featureCollection.features[i] = JSON.parse(result.rows[i].locale);
      featureCollection.properties[i] = JSON.parse(result.rows[i].ttl); //this is wrong
   }

   res.send(featureCollection);
   });

});
}

 function FeatureCollection(){
   this.type = 'FeatureCollection';
   this.features = new Array();
   this.properties = new Object;  //this is wrong
 }


texts:
I am trying to set up GeoServer security to enable external users accessing specific layers from my Geoserver by supplying username and password in eg. QGIS or uDig.

I've connected my "external_user" to an "external_role" and set the
permissions under "Security" -> "data" as:

All_workspaces.*.r for Admin_role

All_workspaces.*.w for Admin_role

External_workspace.*.r for my external_role
Catalog mode: challenge

Under "Authentication", is it possible to use eg. "Form"- authentication and put this above "Anonymous" in the default filter chain? "Form" doesn't seem to be selectable, only Basic which I havn't succeded with so far.
Does anyone have any recommendations for choice of authentication method based on this, basic, form or other?
I've read the docs but I'm new to this and trying to figure out a simple solution that doesn't involve too much advanced configuration.


texts: Where can I find geodata on livestock densities?

texts:
I am currently looking for geodata that represents livestocks (cattle, pig, poultry) at most detailed resolution for France. Up to now, the best data I was able to find is the FAO's GLOBAL LIVESTOCK DENSITIES (census data) data (See FAO data).

For France, the raster has an resolution that is around 3.6km. Does anyone know a better and free to use/not that expensive source? I would be also willing to geocode a livestock farms dataset or something like that.


texts: Return Information from Feature Class Fields by XY Point

texts:
I am trying to get information from fields within specific feature classes based on an XY point.  The user will click on the add-in tool, click on a location in ArcMap, and the tool will find the information from about 5 different feature classes.  This information will be taken and sent to a windows form for the user to see.  I already have code to get the XY click in both screen and map coordinates but don't know how to go about getting the field information.  If I could get 5 strings back to send to the form, that would be perfect. Any help and insight would be great.  Also, feel free to improve my existing code.

Here is what I have so far:

namespace LocationResources
{
    public class LocationIdentify : ESRI.ArcGIS.Desktop.AddIns.Tool
    {
        //Define and initialize global variables
        private string director = "";
        private string city = "";
        private string county = "";
        private string school = "";
        private string subdivision = "";

        private IPoint screenPoint;
        private IPoint mapPoint;

        public LocationIdentify()
        {
        }

        protected override void OnUpdate()
        {
            Enabled = ArcMap.Application != null;
        }

        protected override void OnMouseDown(MouseEventArgs arg)
        {
            if (arg.Button == System.Windows.Forms.MouseButtons.Left)
            {
                try
                {
                    // Get mouse click location in pixels (screen coordinates)
                    screenPoint = new Point();
                    screenPoint.X = arg.X;
                    screenPoint.Y = arg.Y;

                    // Get map document information for transforming pixels to feet coordinates
                    IApplication m_application = Hook as IApplication;
                    IScreenDisplay screenDisplay = (m_application.Document as IMxDocument).ActiveView.ScreenDisplay;

                    // Transform mouse click location into feet (map coordinates)
                    mapPoint = screenDisplay.DisplayTransformation.ToMapPoint(arg.X, arg.Y);
                }
                catch
                {
                    MessageBox.Show("Something went wrong getting the XY point.");
                }

                try
                {
                    // Get spatial information
                    GetSpatialInformation(screenPoint);  // If screenPoint works mapPoint may not be needed
                }
                catch
                {
                    MessageBox.Show("Something went wrong getting the spatial information");
                }

                // Launch location information form
                LocationInformation info = new LocationInformation(director, city, county, school, subdivision);
                info.Show();
            }
        }

        private void GetSpatialInformation(IPoint coordinates)
        {
            // TO-DO --> Get spatial information
        }
    }
}


texts: Efficiency differences between quadtree and geohashes

texts:
I am a newbie to spatial indexing. Can someone please explain me the differences between a quad-tree and geohash. From what I understand both use similar hierarchial methods of dividing the space into 4 quadrants.

However, I am confused as to which one amongst quad-tree and geohash would be efficient for querying and indexing large datasets


texts: Interpreting values produced by GLCM or texture image in ENVI?

texts:
How to use values recorded in stats file ?

Also when i use texture image generated for classification..i am getting error of NaN values.

I have corrected Landsat8 image using Radiometric correction. I didn't perform DOS.

Is this due to this?


texts:
I used QGIS to find the extents/bounding boxes of a list of features from census data.

I used the following method: Vector - Research tools - Polygon From Layer Extent, which creates a new shapefile with a table of extent values for each feature.

However, the resulting extent values (minx, miny, etc.) are not LAT/LONG. I'm not sure what they are. The properties show CRS as "USER:100000 -   Generated CRS (+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0" even though I set the project default to WGS 84 EPSG: 4326. If I change it back to WGS84, the data stays to same.

What CRS format are QGIS extents and how can I convert them to standard Lat/long so that I can plot them?


texts: Split Line by Point (Point-On-Line-Overlayer) for multiple Points

texts:
I am using FME transformer PointOnLineOverlayer to split a line (formed by Road Elements) by Point (Traffic Sign), to do some further processing on split lines. In case there is more than one point on line, my split line is not as expected i.e. for each point it does not create 2 lines, with each line having one end at Point and second end at the end of line.
The attached image should illustrate my problem. Green and red dots represent 2 points (Traffic Signs) on a line (Road). After splitting by FME I am getting 3 lines (double arrowed Black lines) for my line, while I want to get lines as drawn in red and green for each of my points.
Is there any setting in FME transformer to make splitting from point to end of input line?


texts: Is it possible to edit layer transparency in .sld or .qml file for use in QGIS?

texts:
Say I have a shapefile with 20 graduated classes that I'm using in QGIS. I want to save the layer style so that I can use it on other shapefiles, but I can't edit the layer transparency without editing each class (20) individually.

Is it possible to edit the layer transparency in a .sld or .qml file? I arbitrarily set the transparency of one class to 53% and another to 78% (just for the sake of testing), saved both a .sld and a .qml, but I couldn't find these values anywhere in the two files. I also searched for things like 'opacity' and 'transparency' but I don't see an indication of how this is saved. However, my changes (of 53% and 78%) are reflected when loading a .qml file with the shapefiles.

It's not the end of the world to have to edit each class individually at first if the changes are reflected after saving a .qml. But I would still prefer something more reproducible since this would be really easy to do with a simple find/replace, and I may have to change the transparency of all classes at various times in the future.


texts:
I'm sort of newbie to arcpy and trying to create a polyline based on the some predefined points which are put within "network_ND_Junctions1.shp". The code I've written seems like this:

import arcpy

arcpy.env.workspace="C:\\Users\\Erfan\\Desktop\\network_design"
arcpy.env.overwriteOutput="true"

arcpy.CreateFeatureclass_management(arcpy.env.workspace,"Routes.shp","POLYLINE")

route_list=[[2,1,0,4,10],[2,1,0,4,7,16,31,30,28,25,22,10]]

nodes_file="network_ND_Junctions1.shp"

for route in route_list:
    dictionary={}
    for node in route:
        dictionary[node]=[]

    with arcpy.da.SearchCursor(nodes_file,['FID','SHAPE@XY']) as cursor1:
        for row in cursor1:
            if row[0]in dictionary.keys():
                dictionary[row[0]].append(arcpy.Point(row[1][0],row[1][1]))

    myArray = arcpy.Array(dictionary.values())
    myline = arcpy.Polyline(myArray)
    with arcpy.da.InsertCursor("Routes.shp",['SHAPE@']) as cursor2:
        newRow =[myline]
        cursor2.insertRow(newRow)

print myArray.count #output=2
print myline.pointCount #output=0

As you can see in the output comments, arcpy.Array works properly and contains two points but it seems that there's not any point within myline variable. How can I solve this issue?


texts: Can users change symbology of GP output operational layer with Web AppBuilder?

texts:
I built a very simple Python GP script tool in ArcMap 10.4.1 that takes a user input file (e.g. csv) and returns a Point Feature Layer using XY coordinates from the input file.

After I published the results as GP service (no map service returned), I used Web AppBuilder 2.3 (Dev Edition) and set up my GP widget to return the output as an operational layer.

After I run the GP widget and get the result mapped, why am I not allowed/able to change the symbology of my operational layer?


texts: Writing Python script to populate fields of layer with underlying info from another layer?

texts:
I have layers of points and lines that have to be updated regularly.

It needs to have the intersecting attributes from another layer stored in its fields/attribute table.
IE. towns, school districts etc, and roadnames etc.

so for a layer of points, what town, school district and street name each point is in/on.
it was so easy with ArcGIS 9 using VBA, but need a little help getting started using ArcGIS 10.1 and Python.


texts:
I need to find the coordinates for a given location in WGS coordinates. This mean geodetic latitude, longitude, and altitude (with respect to the reference ellipsoid).

I was informed the coordinates, but the altitude was given in MSL (mean sea level).

How do I perform this conversion? Is there a standard way/practice on how to do so?


texts: Moving point labels on vertices inside of polygon

texts:
What I'm trying to do is to label vertices of a polygon with a number and plot that number inside the polygon.

So far I extracted the vertices to a 2nd point layer. Named them correctly and plotted them. I can't find an option to restrict their position to the inside of a polygon. Is there any possibility to do this?

What I know is possible is to create two additional columns in the point layer and use them for plotting coordinates. With lots of points, it can be quite tedious to find the right plotting coordinates for each point.

Is there any cleverer way to accomplish that?


texts: Is there FME transformer to snap point on point of nearest line?

texts:
I have a DWG file which has lines and points. I use FME (2014) to translate it into a shapefile. This works fine but now I want to add some transformers in the workflow.

The first thing I want is to make a simplification of the line network. I did that using Generalizer. Now I want to snap the "points to point on nearest line" of the network.

Is there is transformer for this purpose?

It is a very common tool in GIS packages.

I have read online I could use a combination of "Neighbor Finder" and "2D Point Replacer".

I have followed this workflow:


texts: How to create a built-Up area from multiple polygons?

texts:
I have been struggling over the last few days to emulate an old American map for a city.
I have a shapefile for the buildings of the city, but since they are almost few hundreds, I was trying to combine them together in order to get Built-Up urban areas (blocks) instead. Arcgis Pro provides a feature called Delineate Built-Up Areas but I failed to merge my buildings into different block.
Regards


texts: Calculating a field based on two other fields, if they exist

texts:
I'm using ModelBuilder to iterate through ~800 shapefiles to add and calculate a new field (I'll call it Field_A). Field_A needs to be populated by either Field_B or Field_C, whichever exists. Half of my shapefiles contain Field_B and the other half contain Field_C, but none of the shapefiles contain both.

I tried using the expression: Field_A = !Field_B! or !Field_C! ....but it won't work unless both of those fields are in the attribute table.

How can I get the calculator to recognize whichever field is existing and disregard the other field if it doesn't exist?

(I'm working in ArcMap 10.1)


texts: ogr2ogr does not set osm_id for all multipolygons

texts:
I used ogr2ogr to convert from osm data (pbf) to spatialite :

ogr2ogr -progress -f SQLite -dsco spatialite=yes france.simple.spatialite france-latest.osm.pbf

then, a large part of the imported data misses osm_id

spatialite> select count(*) from multipolygons;
42287946
spatialite> select count(*) from multipolygons where osm_id is null;
42153385

any ideas why, and how to fix this ?


texts: Specify default width and height used by Geoserver's layer preview page

texts:
Currently whenever I try to preview any WMS layer under Geoserver's layer preview page &width=768&height=339 automatically gets set in the request. I'd like to change these values.

Is there a way to specify the default width and height that the Geoserver layer preview page uses?


texts:
This question was just closed for being "opinion-based."  I can't understand how anyone would think this is opinion based.  I'm not asking for the best tracker, or a recommendation for a tracker.

I'm simply asking if anyone knows if a tracker with a particular feature
exists.

That is not a matter of opinion.  It is simply a question of fact.
And what is the feature that I am looking for?  The feature is a tracker that stores its data locally inside the tracker (as opposed to transmitting its data (via cell service or similar) to a cloud storage somewhere).  Does any company manufacture such a tracker?
Again, this is a simple question of fact.  Does such a tracker exist or not.  And if so, which company or companies manufacture it.
I guess some people are quick to close questions, but slow to actually read them!

I am looking to buy a GPS tracking device that I can put in the car that will track (at least) location and speed (perhaps once a minute or so) and store the data.  I do not want to pay for a service to transmit the data in real time.  Instead, I want to just retreive the device periodically and perhaps plug it into my PC (via usb) and pull the data off the device.
I've been searching online for a while and can't seem to find such a device.  Every device I find requires cell service to send the data real time to a server somewhere.  Even if I google for "GPS tracker with no monthly fee" what do I find: trackers where the price includes a year of service!  (I guess an annual fee is the new no monthly fee!)
Is anyone aware of a device that can store the data instead of transmit it?


texts: MapServer: Symbol (v-line) style that automatically follows line direction

texts:
I created this image:

with:

MAP
  SIZE 1000 1000
  IMAGETYPE png24
  EXTENT -30 -10 60 30
  UNITS DD

SYMBOL
  NAME "v-line-shift"
  TYPE vector
  POINTS
    0   0
    5  10
    10  0
  END # POINTS
  ANCHORPOINT 0.5 0
END # SYMBOL

LAYER  # Simple line
  STATUS DEFAULT
  TYPE LINE
  FEATURE
    POINTS
      5 25
      25 20
      45 20
      35 15
    END # Points
  END # Feature
  CLASS
    STYLE
      COLOR 0 0 0
      WIDTH 2.5
    END # STYLE
    STYLE
      SYMBOL "v-line-shift"
      COLOR 0 0 0
      SIZE 8
      INITIALGAP 1
      GAP 10
    END
  END # CLASS
END # LAYER
END # MAP

and I would like to know if these arrows can automatically follow the direction of this changing line (ie. arrow would be aligned with the line, something like e.g., ---->---->---)

I'm using MapServer 6.4.1.

I enclose here an example of the end result.


texts: Deriving coordinates from Application.Display or MxDocument.ActiveView in ArcObjects?

texts:
I've found two possible ways of converting a MouseEventsArg's coordinates to the coordinates of the current map:

The first is taken from a snippet included with ArcObjects:

protected override IPoint OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)
{
    IPoint screenPt = new Point();
    screenPt.PutCoords(arg.X, arg.Y);

    IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
    mxdoc.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint((System.Int32)screenPt.X, (System.Int32)
}

The second was included in some code I've inherited:

protected override IPoint OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)
{
    IPoint pnt;
    IMxApplication app = ArcMap.Application as IMxApplication;
    IPoint pnt = app.Display.DisplayTransformation.ToMapPoint(e.X, e.Y);
}

As far as I can tell, the first is drawing the transformation information from the ActiveView of an MxDocument, while the second takes it from the Display of the ArcMap.Application.

Is there ever a time when these two codes will return Points with different XY coordinates?

Is one way safer than the other?


texts: Leaflet misinterpreting tiles info from an ESRI hosted Tiled Service

texts:
I am working on showing my ESRI tiled layer to a simple Leaflet web viewer.

The tiled layer is hosted from ArcGIS Online.

I tested both the ESRI Leaflet plugin and the original Leaflet's tyleLayer with Leaflet API.

var my_tiled_layer_1 = L.esri.tiledMapLayer({
    url: "{MY_ARCGIS_URL}/MapServer"
})

var my_tiled_layer_2 = L.tileLayer('{MY_ARCGIS_URL}/MapServer/tile/{z}/{y}/{x}');

The initialization of the map worked as intended, and the map tried to add tiles from the correct URL. However when it searched for the {x} and {y} coordinates, something got messed up. Either option gave me a whole bunch of 404 responses. From the Developer Tool, I found the actual tiles it was trying to load did not exist at all. For example Leaflet tried to load

MapServer/tile/8/155/226

However, the actual Level 8 tiles only range from MapServer/tile/8/16028/10942 to MapServer/tile/8/16201/11189

Is there anything I did wrong?


texts: Locating landfall point of cyclone using ArcGIS Desktop?

texts:
On 12.12.2016 Severe Cyclone Storm "VARADHA" crossed Chennai coast.  Nungambakkam [13 Deg 04min 04 sec N/80 Deg 14 Min 44 sec]is 4.8km away from coast.  At 15.33 hrs IST it recorded 974.8hPa pressure and winds were almost calm [0.4knots]

At airport Meenambakkam the pressure is 976.3hPa and it is 10.3 km away from coast.  Varada would have crossed in between these location.

I want to interpolate pressure and wish to find the exact location of crossing with central pressure.

Which interpolation technique do I have to use?


texts: How can I enable panning using Openlayers and a vector overlay with controls?

texts:
I am creating a map using the following code:

map = new OpenLayers.Map({
    div: "map"
});

Then I load in a geojson polygon layer and add two controls

var highlightCtrl = new OpenLayers.Control.SelectFeature(geojson_layer, {
    hover: true,
    highlightOnly: true,
    renderIntent: "temporary",
    eventListeners: {
                beforefeaturehighlighted: unreport,
                featurehighlighted: report,
                featureunhighlighted: unreport
            }
});

selectCtrl = new OpenLayers.Control.SelectFeature(geojson_layer,
    {clickout: true ,
     renderIntent: "select",
     onSelect: selected,
     onUnselect: function(e){console.log('unselected)}
});

map.addControl(highlightCtrl);
map.addControl(selectCtrl);

highlightCtrl.activate();
selectCtrl.activate();

This layer covers the entire map div, and now I can no longer pan the map. Is there a way that I can enable panning when I have a layer with controls that cover the basemap?


texts:
I am still an arcpy newbie, but I think I have a basic understanding of fundamentals including syntax, variables, modules, etc. My question is regarding the usage of parameters while creating a tool.

I am trying to create a tool that will calculate and summarize the area of a feature layer "levels". I have the script working for the calculation and summarization of the feature layer which outputs as a table, but I am not able to get the parameter inputs for the tool working. Syntax checks out (using PyScripter as IDE) but when I try to run the tool in ArcMap "This tool has no parameters" appears. Ideally, I would like the user to input the feature layer "levels" as a parameter and return the summary table, but a message window showing area summary would be even better if thats possible. Please let me know what might be going wrong with the parameters, or how I could change the output to a message window.

import arcpy

class GeoArea(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "GeoArea"
        self.alias = "Tom"

        # List of tool classes associated with this toolbox
        self.tools = [LevelsCalcTool]

class LevelsCalcTool(object):
    def __init__(self):
        self.label = "LevelsCalcTool"
        self.description = "Calculates and summarizes Levels' area in square meters"
        self.canRunInBackground = False

    def getParameterInfo(self):

        param0 = arcpy.Parameter(
            displayName="Input Features",
            name="in_features",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input"
        )

        params = [param0]

        return params

    def execute(self, parameters, messages):
        """The source code of the tool."""
        arcpy.AddField_management(params,"area","DOUBLE","#","#","#","#","NULLABLE","NON_REQUIRED","#")
        arcpy.CalculateField_management(params,"area","!shape.area@squaremeters!","PYTHON_9.3","#")
        arcpy.Statistics_analysis(params, "out_table.dbf", [["area", "SUM"]])

        return


texts:
Possible Duplicate:
  Comparison of WMS and WFS services

I am trying to overlay WFS over WMS and still not successful. So wanted to know when and why would someone use WFS?

I need to extract the info from the WFS layer that I am trying to overlay onto the WMS. If yes how?


texts: Page definition query not turning layer on and off when driven by Definition Querried duplicate layer of index layer?

texts:
I am trying to have some labels turn on and off dynamically with the data driven page being displayed, but it is not working.

First, I have a data driven index layer moving from map to map, called "DD_Index".  This is working.  It is a list of parks driven off the field "PARK"

Second I have a second layer of the FC used for the data driven index, named "Little Park Label Layer", definition queried to only show a single record, in this case a specific park named "Little Park", which is a value in the data driven index field "PARK".  This layer also has a Page Definition enabled with the Page Name Field set to "PARK" and "Match" enabled

Third I have an annotation group in the MXD called "Little Park Labels", it is set with Associated Layer of "Little Park Label Layer"

The theory being when the data driven index gets to "Little Park" the page definition turns on the "Little Park Label Layer", which in turn turns on the associated annotation group "Little Park Labels".

I can test turning the labels on and off with the associated layer manually and it works.  I turn on the "Little Park Label Layer" and the labels in the annotation group display.  All good there.

However if I enable data driven pages and go to the page "Little Park" the "Little Park Label Layer" is not turning on, therefore the labels are not turning on.

Is there an issue with using a definition queried duplicate layer of the index to drive this?


texts: QGIS. Dissolve polygons to create new layer with stats. Not working

texts:
I'm using QGIS 3.2.1 (Bonn)

I'm trying to dissolve and create statistics from a group of polygons, but QGIS is not making it. I do not know if I'm missing something or the software has a bug. These are my steps:

1- Select a group of polygons from a layer

2- Use the "Dissolve" tool of the GDAL add-on.

3- Characteristics to dissolve:

Input layer: "Census ratios", check.

Selected features only, check.

Dissolve field, "Dissolve_" (I have created a column with the same attribute "AAA" for all rows for this purpose), check.

Compute stats for attribute, check.
Attribute selected: Total population (it is an interger64 field), check.

Now take a look at the console:  Dissolve_, SUM(Dissolve_) AS sum, MIN(Dissolve_) AS min, MAX(Dissolve_) AS max, AVG(Dissolve_) AS avg FROM layer_name GROUP BY Dissolve_" -f "ESRI Shapefile"

Is using the Dissolve_ field to calculate the stats. Why? I don't know.

Then go, and the results are crearly wrong. I want to summarize the TOTAL_POPULATION field, not the one used to dissolve the polygons (DISSOLVE_).

Is QGIS failing because of a bug? Or i'm doing something wrong?

Note: I've tried with ArcMap and it goes OK, but i need to work with QGIS.


texts: Offsetting Overlapping Routes in Network Analyst

texts:
I've seen some questions tackle this issue, but havent myself found a working solution. Using ArcMap 10.2.2 Advanced with Network Analyst.

I use the Vehicle Routing Problem to route a couple delivery vehicles around, some resultant routes overlap.

Is there anyway to show all the routes at once? From my picture for instance, the intersection of the yellow/red/green lines, I would like to see all 3 travelling down Valley Blvd, instead of the highest ranked route.


texts:
I am using GDAL 2.01 SWIG in C# to create an imagine file.  When I call Create it pre-allocates the entire file, so if I am writing a 80GB file it I have to wait for it to allocate the 80GB file before I can begin writing to it.

Is there a way to create the 'empty' raster and then fill it as needed?

OSGeo.GDAL.Driver driver = OSGeo.GDAL.Gdal.GetDriverByName("HFA");
OSGeo.GDAL.Dataset dataset = driver.Create("C:\\Test.img", 172032, 98304, 1, OSGeo.GDAL.DataType.GDT_Float32, null);


texts:
I have two tables on CartoDB.

The first one (with name "points_table") has 4 columns.

Country Name (name: country)
Lat (name: lat)
Long (name: long)
the_geom (lat + long)

The second one (with name "links_table") has 2 columns.

Country Name A (name: country1)
Country Name B (name: country2)

The second demonstrate links between two countries. I am new on CartoDB. I found how to visualize the first table, but I still cannot find how I can create curve links between countries from table 2.

Edit

SELECT
  a.cartodb_id AS cartodb_id,
  l.country1 || ' - ' || l.country2 AS link,
  ST_Transform(
  Geometry(
  ST_Segmentize(
  Geography(
  ST_SetSRID(
  ST_MakeLine(a.the_geom, b.the_geom),
  4326)),
  100000)),
  3857) AS the_geom_webmercator
FROM points_table a
JOIN links_table l ON a.country = l.country1
JOIN points_table b ON b.country = l.country2

The SQL Views result contain the cartodb_id column, a column with the two countries and a "the_geom_webmercator" with content as "line". If I hover on it, I have "linestring" title.

However, when I try the MAP view, the website stack on "Loading tiles..." and after a few second, I have a internet connection problem. Any ideas how I can solve it?


texts:
I need the global SOC stocks map at 0-30cm from the new SoilGrids2.0:
https://files.isric.org/soilgrids/latest/data/ocs/

I have no previous experience with VRT files. I have simply downloaded "ocs_0-30cm_mean.vrt" and "ocs_0-30cm_mean.vrt.ovr". Honestly, I don't know if I need one, the other, or both.

In R, I load it with the raster package.

soc <- raster("ocs_0-30cm_mean.vrt.ovr") # t/ha
soc
class      : RasterLayer
dimensions : 14509, 39812, 577632308  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 39812, 0, 14509  (xmin, xmax, ymin, ymax)
crs        : NA
source     : ocs_0-30cm_mean.vrt.ovr
names      : ocs_0.30cm_mean.vrt
values     : -32768, 32767  (min, max)

However, the resolution appears to be 1 degree, when it should be 250m. The extent makes no sense (I guess it should be c(-180, 180, -90, 90), and there is no CRS.

I can plot (soc), though I don't know how to change the resolution to e.g. 0.25.
I have also tried the same with "ocs_0-30cm_mean.vrt" (instead of "ocs_0-30cm_mean.vrt.ovr"), but the file does not plot

plot(soc)
Error in rgdal::getRasterData(con, offset = offs, region.dim = c(1, nc),  :
  Failure during raster IO

Please I would appreciate any feedback to make this raster work, with 0.25 resolution and appropriate extent.