texts: While loop with Select by location (Add to selection) is not stopping

texts:
I have problem with my while loop code. It is not stopping, despite that values used to evaluated the condition to stop the loop are actually met (I can see in in my print statement).

Here is the code:

count_selectionBefore = 0
count_selectionAfter = 1

while (count_selectionBefore != count_selectionAfter):
    count_selectionBefore = arcpy.GetCount_management("inputLayer")
    arcpy.SelectLayerByLocation_management("inputLayer", "BOUNDARY_TOUCHES", "inputLayer","","ADD_TO_SELECTION")
    count_selectionAfter = arcpy.GetCount_management("inputLayer")
    print "count_selectionBefore", count_selectionBefore
    print "count_selectionAfter", count_selectionAfter
    if count_selectionBefore == count_selectionAfter:
        break
        print "loop is over"

print "selection done. Records selected: ", count_selectionBefore

I know that the values of count_selectionBefore and count_selectionAfte are the same from reading them from print statements.

If you can see why my loop isn't stopping, please let me know.

I updated my code using the help below and new code looks like this:

#point to two datasets which are going to be used
firetrails = r'H:\NetworkAnalyst\TopologyVer4\firetrail_mb.gdb\firetrail_step4Export'
roads = r'H:\NetworkAnalyst\TopologyVer4\mb_outputs.gdb\roads_O_noZnoPnoC' #r'H:\NetworkAnalyst\TopologyVer4\mb_outputs.gdb\roads_all'

#make feature layer from firetrails
arcpy.MakeFeatureLayer_management(firetrails, "inputLayer")
#make feature layer from roads
arcpy.MakeFeatureLayer_management(roads, "sourceLayer")
arcpy.SelectLayerByLocation_management(r"inputLayer", "BOUNDARY_TOUCHES", r"sourceLayer", "", "NEW_SELECTION")
countBefore = 0
countAfter = 1
while (countBefore != countAfter):
    resultBefore = arcpy.GetCount_management("inputLayer")
    countBefore = int(resultBefore.getOutput(0))
    arcpy.SelectLayerByLocation_management("inputLayer", "BOUNDARY_TOUCHES", "inputLayer","","ADD_TO_SELECTION")
    resultAfter = arcpy.GetCount_management("inputLayer")
    countAfter = int(resultAfter.getOutput(0))
print "selection done. Records selected: ", countBefore

arcpy.CopyFeatures_management ("inputLayer", r'H:\NetworkAnalyst\TopologyVer4\firetrail_mb.gdb\ft_mb_afterIteration')


texts: How to get data source attributes from PostgreSQL and PostGIS layers

texts: I want to get attributes like schema name and table name from PostgreSQL and PostGIS loaded layers in QGIS 3. I've searched in API docs but just found about set data sources, anyone knows how to get this info?


texts: Different results in OD Cost Matrix analysis and Route analysis

texts:
I am using OD Cost Matrix layer in order to find distances between many origins and destinations points. It works fine, but I want to make sure the results make sense, and since OD Cost Matrix can only produce straight lines (and not true route shapes), I tried to calculate some specific origin-destination pairs using a Network Analyst regular Route layer and compare them to the OD Cost Matrix results.

However, in some of the cases the distances seem to be different between the two methods, even when I use the exact same analysis settings and network restrictions. Sometimes it's only a small gap, but in other cases it's probably a whole different route.

Shouldn't the two tools produce the exact same routes?


texts: Creating a distance matrix for line and/or polygon layers?

texts:
At 2.18.16, the QGIS tool Distance Matrix computes pairwise distances between all points, either within one layer or between two layers.  Here's a screenshot for a single layer.  Note that only points (highlighted) are allowed:

I'm looking for a tool or plugin that extends that Distance Matrix concept to include lines and polygons.  For example, suppose you had a single polygon layer depicting lakes and you wanted to know the distance between each lake and every other lake in that same layer (shoreline to shoreline).  Similarly, consider a situation with two layers, one of them point and the other polygon, wherein you needed to know the pairwise distances between every point and every polygon outline.

I considered the NNJoin tool, because it accepts points, lines and polygons.  Unfortunately it does not meet my need because it computes only the single shortest distance for each feature, whereas I'm interested in all distances.

A review of the Processing Toolbox did not turn up anything useful.

I'm considering a workaround that deconstructs lines and/or polygons to their constituent vertices, converts the vertices to points, applies the Distance Matrix tool to those points, records the pairwise point distances and then joins those to the original line/point layers.

That workaround seems complex, so before I begin, I'm asking for help: Is there an easier way to calculate all pairwise distances between features - point, line or polygon - either within a layer or between layers?


texts: Saving Bing Aerial Images out using MapInfo Pro 16.0

texts:
I'm using MapInfo Pro v16.0 and I'm trying to extract the raster data out for later use.

I've saved my final output already as PDF and JPEG but I'd like to know if there is a way in clipping out the area of interest (using Mapinfo v16) and saving it as a raster data for later use.


texts: QGIS switch Visibility Presets using Python

texts:
I am trying to switch visibility preset using python but I cannot find any relevant method in the documentation . Does anyone knows how I can do it?
The VisibilityPresetCollection class is accessible using the following code QgsProject.instance().visibilityPresetCollection().


texts: "Cross-classification and tabulation" tool is not outputting a table

texts:
In QGIS3 I’m trying to use SAGA’s "Cross-classification and tabulation” tool from the processing toolbox, and based on the documentation (http://www.saga-gis.org/saga_tool_doc/6.3.0/grid_analysis_13.html) and this example I found (https://gis.stackexchange.com/a/264242/117412) I’m expecting it to output both a grid and a table. But every time I’ve tried it (with various inputs) it only outputs a grid. And as the below screenshot shows, it doesn’t even show an output parameter for the table.

So why isn’t the "Cross-classification and tabulation” tool outputting a table?


texts: Select by Location multiple layers using ArcPy

texts:
Using ArcMap 10.3, I would like to Select By Location multiple layers using ArcPy. The error message states line 7 is not correct. I can't find the error though. "Water_Body_Boundary" is my source layer.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
bldgLayer = "Water_Body_Boundary"
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    if layer.name != bldgLayer:
        arcpy.SelectLayerByLocation_management(layer, "WITHIN", bldgLayer)

The error message:

  Runtime error Traceback (most recent call last): File "", line 7, in  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 7359, in SelectLayerByLocation raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000840: The value is not a Feature Layer. ERROR 000840: The value is not a Raster Catalog Layer. ERROR 000840: The value is not a Mosaic Layer. Failed to execute (SelectLayerByLocation)


texts: OpenLayers OSM performance issue on mobile browser

texts:
I face a problem with slowdowns when I try to drag the map on my Samsung Galaxy A5 Chrome browser. Map is rendered by OpenLayers 4.6.5 version. So I've made a Leaflet version with the same source link and... everything works just great, there is no slowdowns and other holds... You can check it by yourself on your mobile.

Here is the OL link: http://178.208.82.80:5938/testMap
Here is the Leaflet link: http://178.208.82.80:5938/testMap2

Here is the OL version code:

<!DOCTYPE html>

    <title>TileHosting</title>

    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/openlayers/4.6.5/ol.css">

© CARTO © OpenMapTiles © OpenStreetMap contributors

Here is the Leaflet version code:

<meta http-equiv="Content-Type" content="text/html; charset=UTF-8" />

<link rel="stylesheet" type="text/css" href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />

Leaflet Example

Here's a map of the countries I've either lived in or travelled through for a month or more.

What is this? Is it means that OL is very unproductive? Or it has to be some special setting?


texts: GPS coordinates coming from vehicle connect to streets coordinates

texts:
I have two tables:

"Vehicle" includes GPS coordinates from the moving vehicles and has around 3.000.000 rows
"streetCoordinates" includes GPS coordinates of the streets, where the vehicles were possibly driving. Its length is around 1.800.000 rows

I want to check, whether the vehicles were on the roads or not. I wrote a query but the thing is, it always fails because it is calculating Cartesian Product - it is checking if there is a match for every point from one table to every point in another table.
Is there any other way to come to a conclusion, which GPS coordinates fall on streets and which do not? Is there a way to do this in QGIS or any other software? Is it possible to this graphically? Up until now, I was using BigQuery and Python for data manipulation and QGIS.


texts:
I have problem with Field Calculator, using ArcGIS 10.2. I tried to write a statement but I don't know what is wrong.

I have two columns:

The first contains a number code. Every code represents some text data. For example: "1" means "trees,etc."

The second column (string) is empty and should contain text data.
Could you help me with Python ?


texts: Am I correctly projecting a geometry using a WKID in ArcObjects?

texts:
I would like to take in any valid WKID for a spatial reference and then project a single geometry to it. The following code seems to be working, but (1.) can someone verify that this is the correct way to do this. (2) Can I just put the ProjectionHint parameter in Project5() to 0/false -- thanks

private IGeometry5 ProjectGeometry(IGeometry5 geom, int outSR)
{
    IGeometryServer geomServer = new GeometryServerClass();
    ISpatialReference spatialReference = geomServer.FindSRByWKID("EPSG", outSR, -1, true, true);
    geom.Project5(spatialReference, 0);
}

http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#/Project5_Method/002m000001rz000000/
"The projectionHint parameter contains information about whether a geometry will cross the coordinate system horizons. If it doesn't, parts of the code that check this and intersect the geometry with the horizons can be omitted. This can dramatically improve performance. Use ISpatialReference3::ProjectionHint to initialize this parameter."

C#
public void Project5 (
    ISpatialReference newSpatialReference,
    int ProjectionHint
);


texts: Filtering points in cloud to find representative set

texts:
I have multiple sets of points that represent 'instances' of GPS traces. When ordered and connected, I have a line that represents each trace. I'm trying to find the most average line through this point cloud.

I've tried several approaches - the current one I'm working on is to put all the points together on the same layer in QGIS, and then run it through the SAGA 'Point filter' tool to remove outliers and select a core set of points. I then re-order them (I store ordering information alongside the lat/long of each point), and I end up with a set of points that pretty well represents the route.

Where I'm struggling is in understanding the QGIS parameters for the 'point filter' tool - I can make the parameters I choose work for some routes, but not others. I can't find any documentation on it and I don't quite understand what they do, which makes it hard to tweak them.

Specifically, what unit is the radius in? What does the attribute do? What do the minimum/maximum number of points do? What are quadrants? How does the tolerance affect the selected criteria? How does the tool choose what the maxima and minima are that it applies to tolerance/percentage to?

I've been changing these parameters myself, and haven't really been able to figure out what they do.

Has anyone used this tool before?


texts: Memory allocation problem when running Random forest regression in python

texts:
I have image with 4 bands . I want to see if I can predict the first band using the 3 other bands.
My problem is that  I cannot fit the model because I get the following error message:

MemoryError: could not allocate 469762048 bytes

I thought that maybe I can break the image into subsets and then to run the algorithm , based on this post
memoryerror in Supervised Random Forest Classification in Python sklearn
but I get stuck after the part of creating "flat_pixels".
in my script, i'm using train_test_split in order to predict the data, so I have X_train,X_test,y_train, _y_test. In the given example, the write create the flat_pixels but seems like it doesn't keep the y values . Also, seems like he doesn't have the "fit" part.
This is what I did:
#load image

img=rasterio.open("myImage.tif")

array=img.read()

#create pandas df
df_all=pd.DataFrame(array.reshape([4,-1]).T)
df_all.columns=['band1','band2','band3','band4']
df_all=df_all.reset_index()
df_all

#create index column in order to keep the original index number
df = df_all.dropna(how='any',axis=0)
df.columns=['index','band1','band2','band3','band4']

#define x,y
X = df.iloc[:, 2:]
#y is band 1
y = df.iloc[:,1:2]

#split the data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

here I have tried to split it into subsets in order to reduce the memory allocationsaccording to the post here memoryerror in Supervised Random Forest Classification in Python sklearn  but it seems to me like i'm losing the y values and then how can I predict it?
#convert to numpy

np_X_train=X_train.to_numpy()
np_X_test=X_test.to_numpy()
np_y_train=y_train.to_numpy()
np_y_test=y_test.to_numpy()

#from the post

flat_pixels = np_X_train.reshape(len(np_X_train),3)
flat_pixels_subsets = np.array_split(flat_pixels, 50)

My end goal: I would like to break my image into subsets in order to be able to fit my model and create prediction.
Any other olution to deal with the memoryError that is not through breaking into subsets is welcome as well.


texts: Getting built environment within administrative areas from OSM

texts:
I have a collection of all administrative boundaries (level 8) in a country. I am interested in playing with the shapes of these towns and cities using Python (shapely, osmnx). But I would like to consider only "inhabited" or "built" parts of the towns. Specifically, I would like to extract the polygon of the grey (but dark and light areas from the administrative polygon (purple dotted line) in overpass.

So far I have tried the buildings API without luck. Is there a specific name for this layer? Once I figure out how to extract this layer I could follow a post like How to Intersect between a vector (administrative areas) and a raster layer (slope)? to obtain the intersection, though I am also happy to hear if someone knows of a slick way to to do this all via Python.


texts: Formula for coordinates (lat, lon), azimuth and distance.

texts:
I would like to know how to find a set of coordinates on a small circle some distance from another point.

For example, lets say A = (39.73, -104.98) #Denver

B = (39.83, -106.06) #point approximately 50 nautical miles away from A along the great circle track to C

C = (40.75, -111.88) #Salt Lake City, approximately 323 nautical miles away from A

How do I get the set of orange coordinates when I don't know the angle but I do know the cross track distance?

For concreteness, lets say the outer points near B have a cross track distance of +- 30 nautical miles and the inner points have a cross track of +- 15nm.

And the points near C have a cross track of +- 25 nm.

Also I have an initial bearing from A to C of -1.34 #radians


texts: count points in polygon returns wrong counts in QGIS 3.0.1

texts:
In QGIS 3.0.1, the vector analysis algorithm 'Count points in polygon' is returning wrong counts for some polygons.  The counts returned for all of the polygons which do not contain any points appear to be correct - i.e. they are all 0.  The counts returned for many, but not all, of the polygons that do contain points are greater than the actual counts (based on my own visual inspection of the points layer rendered on top of the polygons layer) of the points contained in each of the polygons.

I've zoomed in to all of the points contained in the polygon (highlighted in red in the image) to ensure that I count all points, and there are only 16 points.  The count indicated by the 'NUMPOINTS' attribute in the 'Identify Results' pane in the bottom right is 34.


texts: Workaround for deleting elements of shapefiles properly

texts: I realise this is probably a long shot given the amount of discussion around this issue (Deleted polygons seem to revert to older versions,"Phantom" features showing up in QGIS 2.4,edits to attribute table in qgis do not show up in arcgis ) but I wonder if anyone has found a workaround to get Qgis to delete shapefile elements properly (see bug 11007) http://hub.qgis.org/issues/11007 that does not involve saving the shapefile as a new layer. I have over 20 virtual fields attached to the layer in question and so 'save as' would cause me to have to re-write them all each time I edit.


texts: Converting QGIS Processing Modeler Python (.py) file to (.svg)

texts: A Processing Modeler can be saved in QGIS 2.18.15 as a Python (.py) file. Adhering to good advice by @Spacedman in an earlier thread, I wish to know what are the many ways the (.py) files for the QGIS Processing Modeler can be converted into other formats such as the Scalable Vector Graphics (.svg) which can be used in a SVG-editor e.g. Inkscape as suggested by @Spacedman (or for other file formats) where I can create an aesthetically neat graphical representation of the QGIS Processing Modeler, unlike in the QGIS Processing Modeler user interface where it can be challenging to align the elements in the modeller, since there is no Grid (like in yED graph editor).


texts: Drawing irregularly shaped boundary polygon around set of LAS points using ArcMap

texts:
I am working with a set of LAS data which is irregularly shaped, picture a starfish as the point coverage, surrounding the starfish may be water on one side, may be more land elsewhere but no points were collected for it.  In fact most of the areas around this starfish of points is land without points collected, and inside the starfish shape is a little lake.

I want to make a DEM of just the area where points were collected, and there are many ways of doing this, but most will interpolate all the no-point areas inbetween "arms" of the starfish data point pattern.  Making the DEM by using a Terrain Dataset allows for the input of a bounding shape, to theoretically omit no-point areas from being interpolated and only geoprocessing inside of the boundary shape.

However, I have difficulty coming up with a precise bounding polygon of the outside edge of the starfish to use as my geoprocessing limiting extent.

I found a process online in a blog by Clayton Crawford which was picked up by Esri and used as their help file for the process of turning LAS to DEM's, but it involves downloading a custom VBA script and several steps which went a bit over my head.

Does anyone have a suggestion on how to get that starfish pattern of LAS datapoints inside a bounding starfish-shaped polygon, without the use of VBA, or alternatively if there's a way to go right from LAS to DEM raster with no-point areas end up as NoData?

I'm using ArcGIS 10.1 for Desktop.


texts:
I am new to using PostGIS and CartoDB, and I am doing some test SQL inserts into both using PHP.  I've had no problem with points, but polygons are proving to be a bit tricky - my SQL insert works just fine with PostGIS on my local server (or so it seems to me), but it does not work with CartoDB.

Here's my PHP...

<?php
    $cartodb_username = "***";
    $api_key = "***";

    $cartoDBsql = "INSERT INTO polygon_test (poly_geom) VALUES (ST_GeomFromText('POLYGON((-71.1776585052917 42.3902909739571,-71.1776820268866 42.3903701743239,-71.1776063012595 42.3903825660754,-71.1775826583081 42.3903033653531,-71.1776585052917 42.3902909739571))', 4326))";

    $ch = curl_init( "https://".$cartodb_username.".cartodb.com/api/v2/sql" );
    $query = http_build_query(array('q'=>$cartoDBsql,'api_key'=>$api_key));
    curl_setopt($ch, CURLOPT_POST, TRUE);
    curl_setopt($ch, CURLOPT_POSTFIELDS, $query);
    curl_setopt($ch, CURLOPT_RETURNTRANSFER, TRUE);
    $result_not_parsed = curl_exec($ch);
?>

This code does insert points correctly, though, for example using:

$cartoDBsql = "INSERT INTO geometry_test ( the_geom ) VALUES (ST_GeomFromText('POINT($longitude $latitude)', 4326))";

... but polygons do not work.  I imagine that I am not formatting the data in my insert statement correctly, but I am not sure what I am doing wrong.  I've searched the CartoDB documentation and can't find an example of an ST_GeomFromText insert for a polygon.


texts: How to convert point coordinates to point geometry in a different spatial reference coordinate system?

texts:
I have point coordinates returned with wkid 3857 and need to convert to 4326 to get point geometry using PostGIS. I know I can return geometry of the same type using the following:

ST_SetSRID(ST_MakePoint(lat,long),4326)

How can I convert points from a different spatial reference system?


texts: Finding tool to divide polygon into specific areas using ArcPy?

texts:
I am using ArcGIS 10.3.1 and am trying to divide a polygon into several smaller polygons based on a percentage.

I have found a theoretical answer Dividing polygon into specific sizes using ArcGIS?

I am trying to find the tool (code) described by the user whuber or something similar, but all links found on google seem to be dead.

The suggested commercial tool is not an option for me.

Is there a working link to the code I am trying to find or some other tool that would do the same?


texts: How to find all the geometries in a table which are at a particular distance from a given point in postgis?

texts:
I have a table called cartlocations:-

Now I have point geometry of latitude and longitude value, say 23.4877191,86.6927571, I want to find all the geometries which are within a distance of 1000 meters from this given point.
I tried the following query:

select locationid
from cartlocations
where ST_DWithin(ST_GeomFromText('POINT(23.4877191 86.6927571)'),location,1000);

But I was getting this error- HINT:  No function matches the given name and argument types. You might need to add explicit type casts.


texts: How to handle multiple layers for get an unify attribute table?

texts:
I am working with other people making shape files with the same attributes on QGis.
I need an unify attribute table of this group of layers in csv for processing the data.

For now we add the shapefiles to a common folder and we merge the new layers with a unify vector layer. This makes problems when someone are using the unify layer, with simultaneous merges or if in any process the unify layer get broken.
I am new with GIS and I don't know all tools, There is other way for solve this or for create an attribute table of multiple layers?

The solution could be simple with some tool, but my knowledge of QGIS is limited.


texts: Importing feature class into ArcGIS Server feature dataset in geodatabase without stopping service?

texts:
I have a feature dataset which holds data published through ArcGIS Server 10.3.1 and I need to import new datasets into the feature dataset but cannot due to locks.

Do I really need to stop all services which read data from this feature dataset to import a new feature class or is there another way?


texts: Background color changes when projecting a raster with ArcGIS 10.4

texts:
I'm using the Project Raster tool via arcpy:

 arcpy.ProjectRaster_management(
     in_raster = clipped_raster,
     out_raster = projected_clipped_raster,
     out_coor_system = COORDINATE_SYSTEM_WGS_1984_WEB_MERCATOR,
     resampling_type = "CUBIC",
     geographic_transform = TRANSFORM_US_IFR_TO_WGS_1984_WEB_MERCATOR,
     in_coor_system = COORDINATE_SYSTEM_US_IFR)

Before the image has a white background:

But after the projection, its grey:

Any idea why this would be happening?

---- EDIT ----

Per the suggestion by @Fetzer, I check the Stretch setting.  I see it has Percent Clip:

However, when I set it to None per his suggestion, the grey actually gets darker:


texts: How to check if area of interest has clouds in a sentinel scene

texts: I am using python script to download sentinel 2 imagery. given area of interest and cloud condition, it's easy to download the scene. but I am more interested to know if there is some way to check whether clouds are present or not in the area of interest. E.g. even with scene with 90% cloud coverage, there might be a case where the clouds may not be present in the AOI making the scene usable.


texts:
I have a layer containing line segments which are connected by their endpoints (e.g. as within a street network).
I'd like to merge (union) polylines within this layer into one feature until following restrictions are met:

total length of the merged feature is greater than a minimal length (e.g. 100m)
only line segments can be merged which are neighbors

The length of the lines are defined within attributes.
It shall be possible to remember which lines are merged (mapping to the original feature ids).

The idea is to iterate through all features with length < minlength and intersect with other features to find the neigbours. Then the neighbour with minimal length will be merged (union) with this feature.
The importat part would be not to iterate through features already merged (maybe by recursion).

Pseudo Code:

Program Mergelength
define minLength
define mergedFeatureList = empty
iterate through all features of featureClass where length < minLength
   if feature not in mergedFeatureList
      select features from featureclass which intersect (touch) feature
      find feature with minimum length
      union feature geometries
      update feature length
      insert feature into mergedFeatureList
      do (repeat) until feature length > minLength

Incomplete code (qgis python):

from PyQt4.QtCore import *
import math
##################################################################################################
minLength = 50

inlayer = qgis.utils.iface.activeLayer()
lenidx = inlayer.fieldNameIndex('REALLENGTH')
lididx = inlayer.fieldNameIndex('RSN')

# Function definition is here
def intersection(feature, layer, idIdx, lenIdx, alreadyMergedSegements):
    intersections = {}
    for ifeature in layer.getFeatures():
        if feature.geometry().intersects(ifeature.geometry()) and ifeature.attributes()[idIdx] not in alreadyMergedSegements:
            intersections[ifeature.attributes()[idIdx]] = ifeature.attributes()[lenIdx]
    return intersections

# Function definition is here
def mergeSegements(lineFeature, lenidx, lididx, minLength, inLayer):
    alreadyMergedSegements = []
    myIntersections = intersection(lineFeature, inLayer, lididx, lenidx, alreadyMergedSegements)
    minRsn = min(myIntersections, key=lambda k: myIntersections[k])
    request = QgsFeatureRequest().setFilterExpression( u'"RSN" = ' + str(minRsn))
    for feature in inLayer.getFeatures(request):
        fet = QgsFeature()
        unionGeom = lineFeature.geometry().combine(feature.geometry())
        mergedLength = float(feature.attributes()[lenidx]) + float(lineFeature.attributes()[lenidx])
        alreadyMergedSegements.append(minRsn)

    return -1

print lenidx
iter = inlayer.getFeatures()
for feature in iter:
    length =  feature.attributes()[lenidx]
    if length < minLength:
        feat = mergeSegements(feature, lenidx, lididx, minLength, inlayer)
    else:
        print "no merge needed."


texts: error: Gdalbuildvrt does not support rotated geo transforms?

texts:
I try to convert several tiff to vrt -file. All tif-files are in same projection (EPSG:3067). I get error message: Gdalbuildvrt does not support rotated geo transforms. What this mean and how to avoid this error? I suppose there isnt any rotation in those tif-files.

Tiffinfo gives next information about example file:

Z:\A_tuotteet\Osoitekartta_40k\ETRS_TM35FIN>tiffinfo Kuvaruutu\Tiff_KESA\K241.tif

TIFFReadDirectory: Warning, Unknown field with tag 33550 (0x830e) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 33922 (0x8482) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34735 (0x87af) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34736 (0x87b0) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34737 (0x87b1) encountered.
TIFF Directory at offset 0xbcafa (772858)
  Image Width: 12000 Image Length: 6000
  Resolution: 254, 254 pixels/inch
  Bits/Sample: 8
  Compression Scheme: LZW
  Photometric Interpretation: palette color (RGB from colormap)
  FillOrder: msb-to-lsb
  Orientation: row 0 top, col 0 lhs
  Samples/Pixel: 1
  Rows/Strip: 2
  Planar Configuration: single image plane
  Page Number: 0-1
  Color Map: (present)
  White Point: 0.3127-0.329
  PrimaryChromaticities: 0.640000,0.330000,0.300000,0.600000,0.150000,0.060000
  Tag 33550: 4.000000,4.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,20000.000000,6642000.000000,0.000000
  Tag 34735: 1,1,0,8,1024,0,1,1,1025,0,1,1,1026,34737,22,0,2049,34737,7,22,2054,0,1,9102,2062,34736,3,0,3072,0,1,3067,3076,0,1,9001
  Tag 34736: 0.000000,0.000000,0.000000
  Tag 34737: ETRS89 / TM35FIN(E,N)|ETRS89|

Same information using gdalinfo

Driver: GTiff/GeoTIFF
Files: Kuvaruutu\Tiff_KESA\K241.tif
Size is 12000, 6000
Coordinate System is:
PROJCS["ETRS89 / TM35FIN(E,N)",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",27],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3067"]]
Origin = (20000.000000000000000,6642000.000000000931323)
Pixel Size = (4.000000000000002,-4.000000000000001)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=254
  TIFFTAG_YRESOLUTION=254
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   20000.000, 6642000.000) ( 18d28'18.40"E, 59d38'16.57"N)
Lower Left  (   20000.000, 6618000.000) ( 18d31'33.02"E, 59d25'29.24"N)
Upper Right (   68000.000, 6642000.000) ( 19d18'54.73"E, 59d41'25.43"N)
Lower Right (   68000.000, 6618000.000) ( 19d21'50.59"E, 59d28'36.50"N)
Center      (   44000.000, 6630000.000) ( 18d55' 8.15"E, 59d33'29.37"N)
Band 1 Block=12000x2 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)


texts: error using gdal_translate with -projwin option

texts:
I have to crop geotiff image, I'm using gdal assembly and I'm trying to use the gdal_translate method

If I crop my image using the -srcwin option all is fine

gdal_translate -srcwin 0 0 400 400 imageBase.tiff image_cropped.tiff

I have to manage bbox in wgs84 coordinates too, so I try to use -projwin option:

gdal_translate -projwin 103.26 23.02 103.36 22.9 imageBase.tiff image_cropped.tiff

I get the following error:

Input file size is 450, 1047
The -projwin option was used, but the geotransform is
rotated.  This configuration is not supported.

My geotiff image has been created programmatically in C# from Jpegcovnerter using tiepoints calculated with my feature image fooprint.

This is the gdalinfo over my geotiff file

Size is 450, 1047
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
GeoTransform =
103.2667570026355, 0.0005062068032082341, -1.832943060114165e-007
23.02489610024383, -5.790469882757381e-005, -0.0005021561705462942
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=JPEG
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 103.2667570,  23.0248961) (103d16'0.33"E, 23d 1'29.63"N)
Lower Left  ( 103.2665651,  22.4991386) (103d15'59.63"E, 22d29'56.90"N)
Upper Right ( 103.4945501,  22.9988390) (103d29'40.38"E, 22d59'55.82"N)
Lower Right ( 103.4943582,  22.4730815) (103d29'39.69"E, 22d28'23.09"N)
Center      ( 103.3805576,  22.7489888) (103d22'50.01"E, 22d44'56.36"N)
Band 1 Block=450x1047 Type=Byte, ColorInterp=Red
Band 2 Block=450x1047 Type=Byte, ColorInterp=Green
Band 3 Block=450x1047 Type=Byte, ColorInterp=Blue

I didn't implemented the class building the Tiff so it looks like my geotiff doesn't support the gdal_translate -projwin option.

Do someone have a idea about this problem and how to solve it?


texts:
I am trying to analyze the MODIS EVI 250m, 16day dataset. I found the R MODIS package and am following the instructions here: https://conservationecology.wordpress.com/2014/08/11/bulk-downloading-and-analysing-modis-data-in-r/

But when I try to run the runGDAL function it tells me that:

library(MODIS)
dates1 <- as.POSIXct( as.Date(c("01/01/2015","31/12/2015"),format = "%d/%m/%Y") )
dates2 <- transDate(dates1[1],dates1[2])
product <- "MOD13Q1"
h = "25"
v = "07"
runGdal(product=product,begin=dates2$beginDOY,end = dates2$endDOY,tileH = h,tileV = v,SDSstring = "111", outProj="4326")

However I get the error:

  Error in runGdal(product = product, begin = dates2$beginDOY, end = dates2$endDOY,  :
  GDAL not installed or configured, read in '?MODISoptions' for help

I tried using the checkTools() command to set GDAL:

MODIS:::checkTools("GDAL")
Checking availabillity of GDAL: OK, GDAL 1.11.3, released 2015/09/16 found!

Does anybody have any ideas on how to fix this issue?

I'm using MODIS version 0.10.33. Output of MODISoptions():

STORAGE:
_______________
localArcPath : /Users/sujoybhattacharyya/MODIS_download
outDirPath   : /Users/sujoybhattacharyya/MODIS_download/Processed
DOWNLOAD:
_______________
MODISserverOrder : LPDAAC, LAADS
dlmethod         : auto
stubbornness     : high

PROCESSING:
_______________
GDAL           : Not available. Use 'MODIS:::checkTools('GDAL')' for more  information!
MRT            : Not available. Use 'MODIS:::checkTools('MRT')' for more information!
pixelSize      : asIn
outProj        : asIn
resamplingType : NN
dataFormat     : GTiff

DEPENDENCIES:
_______________


texts: Finding number of left and right turns for each shortest route

texts:
I have a network dataset for a city. I have, lets say, 200 trips that are made in that city. I constructed the shortest route for the every pair of the trips, and also the actual route travelled by them using the GPS points available.
I want to include the number of turns, either left or right, for the shortest route as well as observed route for all of the 200 trips.
Is there any tools in network analysis to do this for a large dataset, that doesn't requires a huge labor and time? (I have thousands of trips data so I was looking for a simple tool or function to accomplish this)

I tried using the COPY TRAVERSED SOURCE FEATURES tool, as I found that this could give the required result for my purpose. However, the output for the turn feature is empty every time I run this tool. Why the table for the TURN feature is empty? What might be the reason.


texts: Adding third party online basemap in ArcMap

texts:
I am trying to figure out the coastline changes between past and present using map of a Japanese web page, something like Google Earth.

Is there any way I can import the basemap of an online web page into ArcMap so that I can trace the coastline to compare if there are any changes between past and present?


texts: How to make transparent WMS vector background in QGIS?

texts:
I have a WMS vector layer added to QGIS desktop.
But the background is not transparent - it is white.
How to change that?
When I add this WMS to ArcGIS - all is just fine.

Here is the WMS link:

http://maps.rosreestr.ru/arcgis/services/Cadastre/CadastreWMS/MapServer/WMSServer?


texts: Split tool not splitting line as expected in ArcMap

texts:
I have a long line that I would like to split. Therefore I start editing, select the line and choose the split tool. I then click where I want to divide the line. What I notice is that on some occasions what I get is not what I would like, a line divided into two segments, but I get a line split into three segments where 2 of them are connected to each other (as in the image below)

How can I solve this problem?


texts: Speed differences between SQL Server 2012 and PostGIS UnionAggregate

texts:
There seems to be a massive difference in speed between PostGIS and SQL Server 2012 with a union aggregate function.

SELECT geometry::UnionAggregate(the_geom)
FROM loc
GROUP BY LocalityCo

vs

SELECT ST_Multi(ST_Union(the_geom))
FROM dummydata.loc
GROUP BY "LocalityCo"

I get less then a second in SQL Server 2012 but around 25 seconds in PostGIS.

PostGIS version is: 1.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

Is there anything I can do to speed up the PostGIS query or it's just the way it is?

Note: Both have spatial indexes.


texts: Web AppBuilder for ArcGIS Enterprise login popup for each refresh

texts:
I have created a Web AppBuilder web application and downloaded. And published it on IIS. When I visited the application via browser the login popup is appearing like following.

But this popup is appearing every time when I refreshed the page. The sign in info does not saved on browser.  If I do Map share settings for everyone, the popup does not fire. I actually hope to redirect me portal login page. How can I solve this problem?


texts: Performance problems using big points shapefile in TileMill

texts:
I've downloaded all the addresses in Italy from OSM and now I've a shapefile with more 514.000 points (16.065 KBbyte when zipped).

I need to publish them with other geospatial information (points ...): I've not a public GIS infrastructure with some GIS engine (GeoServer, MapServer, etc.), so I've planned to use TileMill to create tiles and then to use a little Leaflet page to publish all togheter.

At the moment I'm not able to generate tiles because TileMill is very very slow to upload the shapefile (I've started with OSM addresses that is the biggest one ....), and it's quite impossibile to export them in mbtiles. Note that the shapefile has the .shx and is shown right in QGIS. I think that I'm wrong somewhere Any suggestion? Alternatives?

I'd to use both raster tiles and vector tiles when I'll show the details.

You can see my points data at this address https://www.dropbox.com/s/a32pudxht18z832/CiviciOSM.zip


texts: What does st_hausdorffdistance actually compare?

texts:
Does st_hausdorffdistance just look at discreet points used to construct a line segment or does it compare two line segments as a continuous set of points?

Also, the points used to construct my line segments are in degrees so the result of st_hausdorffdistance returns a number in degrees as well which is not really sensical.

How can I convert that to meters?


texts: Cached Featured in WFS - can I change default

texts: QGIS currently sets the default to cached download for adding WFS - is it possible anywhere to set the default to non cached mapping (I know you can untick the box, but its for a larger user group)?


texts:
I created a model which geocodes addresses using an excel file. I wanted to know if there is a way I can give a different excel file every time to the model. I do not want to manually add a new excel file every time.

Can a python script point to a folder with multiple excel files and give them, one by one as an input to ModelBuilder?


texts: Maximum number of shapefiles for merge tool?

texts:
I have over 180,000 shapefiles that I have to merge (or append) into a final shapefile.
I have written a python script to do so which uses the acrpy.append_management tool. However the runtime is too long.

I was thinking of maybe testing out the merge tool but I am skeptical that it will accept this many shapefiles and merge them together in a reasonable amount of time.

Does anyone have any experience in merging large datasets?
Is there an efficient way to merge large amounts of shapefiles?


texts:
I am Looking for a python script that will automatically assign a unique sequential value to a field every time a new feature (asset) is created. This ID will also need the abbreviated name prepended to the identifier.

So for example, if I add a new sign (SN) to the signs feature class, I'd like that sign to have an Asset ID of SN235.

Alternatively, I could field calculate each time a new asset or series of assets is added so long as the next feature is given the next sequential value.

I have seen scripts such as this: Auto-incrementing in Python Script with ArcPy cursor?

But it isn't quite what I need.

Has anyone created or come across a script that will do this?


texts:
Kind of a newbie with ArcObjects .NET. Trying to open a dockable window with the following code:

public class CAD_Test_01_button : ESRI.ArcGIS.Desktop.AddIns.Button
{
    public CAD_Test_01_button()
    {
    }

    protected override void OnClick()
    {

        UID dockWinID = new UIDClass();
        dockWinID.Value = ThisAddIn.IDs.DockableWindow1;

        // Use GetDockableWindow directly as we want the client IDockableWindow not the internal class
        IDockableWindow dockWindow = ArcMap.DockableWindowManager.GetDockableWindow(dockWinID);
        dockWindow.Show(true);

    }

    protected override void OnUpdate()
    {
    }
}

But I get the errors:

Error: The type or namespace name 'IDockableWindow' could not be found (are you missing a using directive or an assembly reference?)

Error: The name 'ArcMap' does not exist in the current context

What references could I be missing?

Thx,

Rich


texts:
I'm using jp2 with no problem for reading and can create a copy with createcopy, but when I try to create jp2, it fails.

driver = gdal.GetDriverByName("JP2OpenJPEG")
dst_ds = driver.Create("test.jp2",xsize=img.RasterXSize, ysize=img.RasterYSize,bands=1,eType=gdal.GDT_UInt16)

Same code but using GTiff driver works ok:

driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create("test.jp2",xsize=img.RasterXSize, ysize=img.RasterYSize,bands=1,eType=gdal.GDT_UInt16)

In the first case dst_ds is none, in the second one I get a correctly created dataset.

Note: img is another JP2 file opened without problem for reading.

Originally I was doing a createcopy, update some pixel values, doing a  WriteArray, and flushing but I can see the changes reflected in de destination file, so I decide to create from scratch.

Is there a write limitation in the driver? I'm doing something wrong?

UPDATE 1:

adding the code:

if(dst_ds is None):
   print(gdal.GetLastErrorMsg())

I'm getting the message:
GDALDriver::Create() ... no create method implemented for this format.

It does not look good....


texts: How to implement close Popup for Multipolygon?

texts: I am closing the popup for geoJson (layer.closePopup()) to make a modal to appear. The function works fine for regular geometry but fails for Multipolygon. I get the error - ClosePopup is not a defined function.


texts: Command-line utility to convert GPX into KML?

texts:
I use Cycle.travel to draw cycling routes. However, the Android application I use on my smartphone only supports KML.

As an alternative to GPSVisualizer, I tried the command-line utility GPSBabel, but, besides turning a compact 5K GPX into a much bigger 80K KML, the route is displayed as a long series of waypoints instead of a simple line:

Does someone know of an alternative command-line utility to GPSBabel?

Here's a sample from the source GPX file and the KML generated by GPSBabel:

GPX

<trk>
<name>ACTIVE LOG</name>
<trkseg>
<trkpt lat="48.81569" lon="2.31610" />
<trkpt lat="48.81579" lon="2.31590" />
…
<trkpt lat="48.75456" lon="2.18915" />
<trkpt lat="48.75510" lon="2.18814" />
</trkseg>
</trk>

KML from GpsBabel

…
    <LookAt>
      <longitude>2.252120</longitude>
      <latitude>48.785320</latitude>
      <range>15131.199708</range>
    </LookAt>

…
…
<![CDATA[

Distance 9.8 mi

]]>
</description>
        <Folder>
          <name>Points</name>
          <Placemark>
            <name>ACTIVE LOG-0</name>
            <snippet/>
            <description><![CDATA[

Longitude: 2.316100

Latitude: 48.815690

Heading: 2.0

…
        <Placemark>
          <name>Path</name>

          <LineString>
            <tessellate>1</tessellate>
            <coordinates>
2.316100,48.815690
2.315900,48.815790
…
2.189150,48.754560
2.188140,48.755100
</coordinates>
          </LineString>
        </Placemark>
      </Folder>
    </Folder>
  </Document>
</kml>


texts: How to determine the lat long of a point in the atmosphere which is 250 km altitude?

texts:
Regarding determination of coordinates of Ionospheric Pierce Point (IPP)

S is the position of a geostationary satellite. E is the elevation angle.

The radio signal coming from the satellite s penetrates the IPP point and received at location whose lat long are 16.8 N and 74.2 E.

We need to find the coordinates (lat, long) of a point (IPP) which is in the atmosphere.

We need the latitude and longitude of that point, which is called the ionospheric pierce point (IPP).


texts: Is there a difference between the Planchon & Darboux and Wang and Liu depression filling algorithms? Other than speed?

texts:
Can anyone tell me, based on actual analytical experience, if there is a difference between these two depression filling algorithms, other than the speed at which they process and fill depressions (sinks) in DEMs?

A fast, simple and versatile algorithm to fill the depressions of digital elevation models

Olivier Planchon, Frederic Darboux

and

An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling

Wang and Liu

Thanks.


texts:
I have a QTreeWidget in a QWidget which is the gui of my plugin(cpp). Now i want to get the keyPress-Event of a Node in the TreeView. I try to reimplement the signal of the QWidget.

protected:
    virtual void keyPressEvent(QKeyEvent *);

But this don't work. I think the problem is: My gui class derives from QWidget(the signal will be reimplemented), but it has also added the QTreeWidget as member, and the signal from it should be reimplemented? But how to do this? Any suggestions?


texts: Where can I find some Sample Data for Linear Referencing?

texts:
I'm teaching GIS at a Local college, and one of the classes is on Linear Referencing, specifically of Roads.

For the Practical classes, I need some Sample Data for Linear Referencing. Ideally, this would be:

In Shapefile Format
Would be Roads or Highways.
Each Feature would have the staring & ending Mile Measurements.

I have not been able to find any freely accessible data which can be used for such purposes.

Where can I find this kind of data?


texts: OpenLayers Map with custom SRS (without EPSG code)

texts:
Is it possible to show geodata by OpenLayers in custom projection?

I want to show my geodata in OpenLayers interface.

Source layer(s):  shapefile(s), standard epsg:4326 lat/lon projection; processed by MapServer (WMS Layers).

Result map: Albers Conic Equal Area projection, with specific central meridian (it has no EPSG code).

Definition in PROJ.4 format:

+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=93 +x_0=16500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs

As to MapServer Layers – they are all well, I can show geodata through standard MapServer CGI template-based interface, in this custom Albers projection…

For OpenLayers I tried to define Proj4js.defs –

<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<link rel="stylesheet" href="ol/style.css" type="text/css" />

OpenLayers Test

It seems that the presented coordinate transformation works correctly (displaying the appropriate numbers when I move the mouse), but the map does not appear…

How can I do this?

Discussion of the problem & Answer to the question --
It looks like I found the answer myself...  The above source code has been corrected, and already has a working version

I changed WMS-layer upon layer of MapServer type

Original (non-working) version was:

...
    var rus_general = new OpenLayers.Layer.WMS( 'Russian Regions',
        'http://localhost/cgi-bin/mapserv.exe?map=C:/OSGeo4W/apache/htdocs/test/z-albers.map',
    {layers: 'rus_general'}
        );
...

In this case, the bottleneck - WMS-service: I do not know how to specify a custom (non-EPSG) projection in this case...

This is not a problem for the MapServer layer type. MapServer generates a map projection with PROJECTION section in the MAP-file. There are no restrictions such as "EPSG Only"...


texts: Project Tool without an exisiting Geographic transformation using ArcGIS Pro

texts:
I'm trying to insert a layer and in my project and while several other layers have been projected correctly on-the-fly by ArcGIS Pro, the final polygon layer is creating trouble. It slightly misfits to the other layers (approx. 50m away from the proper place, in a 1:5.000.000 scale). I'm aware that this is linked to a projection problem, but I cannot figure out how to solve it.
What I tried so far:

Use the Project Tool with the CRS of the base maps and others
Use the Define Projection Tool with the CRS of the base maps and others
Use 1) and 2) consequently
Trying to create a Geographic Transformation-Oneliner manually, because the drop-down bar never appeared. Maybe the CRS is too exotic or wrongly named?
Creating a placeholder Featureclass with Unknown CRS to first project the Unknown CRS on the problematic layer and then redo 1)
Export Features from the problematic  layer and using the in-tool projection


texts: Shapefile to PostGIS using FME - Column has M dimension, geometry does not

texts:
I'm trying to insert a shapefile into PostGIS using FME but I always get an error saying:
ERROR: Column has M dimention but geometry does not
CONTEXT: PL/pgSQL function debug.trig_fonc_force_geom_3d() line 6 at assignment

If I inspect the features in FME Data Inspector, the geometry info for the existing data in postGIS looks the same as the shapefile data, i.e.:
IFMELine (2 Coordinates)
Closed           No
Measures (1)     <default_measure>
Coordinates (2)  Coordinate Dimension: 3
0                -3.432432, 40.23948723, 0, <0>
1                -3.238473, 40.42348234, 0, <0>

I have tried every combination of setting/removing Z and M values but it doesn't make any difference.  I've tried both Spatial Type options in the wirter (Geometry/Geography) - what can I do?
Note:
I'm inserting into an existing table containing geometry.  When the table was empty I was able to use FME to insert one batch of data but now I'm getting this error.


texts: How to tell if a line intersects a polygon not only at endpoints?

texts:
How do I do this in .net topology suite or jts? Please see image.

Any other library that can I run on .NET is fine, or even a general algorithm if there is no option.

For example this answer is only for endpoints:
Check if line crosses a polygon

Right now I'm randomly experimenting and if i intersect and check

if (intersection.Boundary.ToString().Contains("EMPTY"))

it seems to work, but... yeah no idea why or if it's all cases.


texts: convert headerless flat binary files to tif

texts:
I have set of headerless flat binary files (those are out put files of TIMESAT software). I want to convert those to .tif format.
I did it using ENVI. but it does not overlap with the original input files of TIMESAT(but both input and out put files are shows the coordinate system and datum as WGs84.). in the processing of TIMESAT no coortination transformation is occured.
can any one tel me what is the problem of ENVI.
or do you know some method to convert headerles flat binary to tif.


texts: Filter legend entry for raster layer QGIS

texts:
Is it possible to use the "Filter legend by Map content" button, or achieve the same functionality with a raster layer?

I've got a map that has a dozen different raster layers, but am only displaying one at time in a map composer.  The item properties for legends has a button for filtering legend content to only include legend items for layers that are displayed in the current composer.  This works very well for vector layers, but doesn't seem to have an effect on raster layers.

I would like to include layer entries for each of the raster layers, but only display them when the associated layer is displayed.


texts: Dynamic Contour generation and project on google maps

texts:
Is it possible to generate contour map layer at runtime and overlap or project on google maps.

I am using google map api v3 and want to generate contour map layer using data from sql server and project on google map.


texts: Generating vector contour map of fictional terrain

texts:
Originally I had asked this on the Graphics Design SE but after some realization of the question I was asking, it's probably better to ask it here.

I'm making a fantasy world map in illustrator that I'd like to take to the next level by adding terrain contours. Currently, it's a flat political map like this map of Canada for example. I'd like to add an effect like what was used in this map here to show the relative depths of the water:

This is obviously not done by hand. They would have used a program that compiled real terrain data into vector contours which would get exported to Illustrator/Indesign/editor of choice.

I don't have the luxury of real data, so how would I go about generating my own using the map I've already drawn as a starting point?

I have a good idea of what it is I need to do -- I'd like to export my existing map to a tool that will generate/edit terrain and then export it back as a vector for Illustrator. I'm just totally new at GIS and don't have any idea of the software or tools of the trade. I've looked glancingly at things like R and Artificial Terrain Generation, and have used tools like QGIS in the past, but I could really use some direction on what's the best way to go about doing this and what software programs or libraries could help me.

I'm no stranger to coding and know C and Python pretty well.


texts: OpenStreetMap - Determining number of nodes, ways, and relations in a PBF file?

texts: I'm currently importing a subsection of North America into PostGIS to render my own basemap tiles using osm2pgsql. Is there any way to get an accurate count of the number of nodes, ways, and relations contained within my PBF file that I am generating to determine how long this will take to process on my computer?


texts:
I am seeing mixed answers on other posts so I thought I would ask here.  Is there a way to automate the raster calculator?

I am trying to create a series of inundation maps at fine increments.  It is tedious to do manually, so I would like to write a script that will loop the raster calculator.


texts: Auto-refreshing triggered by changes in directory for ArcPy?

texts:
I have a directory change monitor python script which monitors any changes happened in a directory, I want to trigger Arcmap auto refreshing whenever a change has been detected (eg. a new line added to a file). I tried to launch the monitor program within arcmap python command window using execfile("myprogram.py"). The problem I encountered is that once the script starts to monitor and auto-refresh, the arcmap main window was locked and frozen, no way to stop it except shutting down the arcmap, although I can see the auto-refresh happened whenever there is a change. Is there a way to run the monitor script in a separate process or background without locking or freezing the main window? Here is my monitor script which is borrowed from Tim Golden,

    import os
    import win32file
    import win32con
    import arcpy

    # Search for the new data
    ACTION = {
       1: "Created",
       2: "Deleted",
       3: "Updated",
       4: "Renamed from something",
       5: "Renamed to something"
    }

    FILE_LIST_DIRECTORY = 0X0001
    path_to_watch = "C:\Directory"

    hDir = win32file.CreateFile(
          path_to_watch,
         FILE_LIST_DIRECTORY,
         win32con.FILE_SHARE_READ | win32con.FILE_SHARE_WRITE |   win32con.FILE_SHARE_DELETE,
         None,
          win32con.OPEN_EXISTING,
          win32con.FILE_FLAG_BACKUP_SEMANTICS,
          None
     )

     while 1:
        results = win32file.ReadDirectoryChangesW(
            hDir,
            1024,
            True,
            win32con.FILE_NOTIFY_CHANGE_FILE_NAME |
            win32con.FILE_NOTIFY_CHANGE_DIR_NAME |
            win32con.FILE_NOTIFY_CHANGE_ATTRIBUTES |
            win32con.FILE_NOTIFY_CHANGE_SIZE |
            win32con.FILE_NOTIFY_CHANGE_LAST_WRITE |
            win32con.FILE_NOTIFY_CHANGE_SECURITY,
            None,
            None
         )

    for action, file in results:
        full_filename = os.path.join(path_to_watch, file)
        print full_filename, ACTION.get(action, "Unknown")
        arcpy.RefreshActiveView()


texts: How to create isometric 3D map with postgis

texts:
I am looking for a way to create isometric 3D maps based on data stored in a PostGIS database.

My idea is to create a PostGIS function that will return a SVG image like the one below. It is easy to do for a 2D image using ST_AsSVG and some math to position the map, but I don't know how to "rotate" the map to add some perspective, then "extrude" the contours to add volume. Any help is welcome !


texts: Count points in polygon (weighted) for multiple attributes in QGIS 3 in one go

texts:
Is there a way to count weighted points inside polygons for multiple attributes at one time?

I can count weight points (below) for an individual attribute, which creates a new polygon layer with an updated attribute table. I have ten attributes of the points that I need to aggregate into the polygons, so this would take a long time to do one by one.

I'd considered setting them up as batch processes, but this generates a new layer each time. I'd like to only generate a new layer once.


texts: Troubling using MODIS package in R (mac). Error with gdal

texts:
I'm new to using MODIS data, and have been trying to set up a a way to use the HDF files produced by MODIS on R. I had trouble with gdal (could only download the package through system dialog), and MODIS (had to install directly from tar.gz saved to desktop) packages, but now have run into a problem I can't quite get past. I've been trying to use the MODIS package in R, as directed by the tutorial: https://conservationecology.wordpress.com/2014/08/11/bulk-downloading-and-analysing-modis-data-in-r/

The code so far is as follows:

library(MODIS)
library(rgdal)
dates <- as.POSIXct( as.Date(c("01/05/2014", "31/05/2014"), format = "%d/%m/%Y"))
dates2 <- transDate(dates[1], dates[2])
h = "21"
v = "09"
product <- "MOD13Q1"
runGdal(product = product, begin = dates2$beginDOY, end = dates2$endDOY, tileH = h, tileV = v,)

on running the last line I get the error

########################
outProj          =  asIn
pixelSize        =  asIn
resamplingType   =  near
Output directory =  /Users/Mapes/MODIS_ARC/PROCESSED/MOD13Q1.006_20160403221952  (no 'job' name specified, generated (date/time based))
########################
sh: gdalinfo: command not found
Error in system(paste("gdalinfo ", HdfName, sep = ""), intern = TRUE) :
error in running command

so that when trying the next command:

vi <- preStack(path="/Users/Mapes/MODIS_ARC/PROCESSED/MOD13Q1.005_20140810192530/", pattern="*_EVI.tif$")

the error:

No files found!

appears.

Checking MODIS:::checkTools("GDAL") only shows

Checking availabillity of GDAL:
OK, GDAL 1.11.4, released 2016/01/25 found!

and MODSOptions() outputs:

To install all required and suggested packages run:
setRepositories() # activate CRAN, R-forge, and Omegahat and then:
install.packages(c(' ptw '),dependencies=TRUE)

'MRT_HOME' not set/found! MRT is NOT enabled! See: 'https://lpdaac.usgs.gov/tools/modis_reprojection_tool'

STORAGE:
_______________
localArcPath : /Users/Mapes/MODIS_ARC
outDirPath   : /Users/Mapes/MODIS_ARC/PROCESSED

DOWNLOAD:
_______________
MODISserverOrder : LPDAAC, LAADS
dlmethod         : auto
stubbornness     : high

PROCESSING:
_______________
GDAL           : GDAL 1.11.4, released 2016/01/25
MRT            : Not available. Use 'MODIS:::checkTools('MRT')' for more information!
pixelSize      : asIn
outProj        : asIn
resamplingType : NN
dataFormat     : GTiff

DEPENDENCIES:
_______________

>

Is there something I'm missing? I'm still not quite sure why the runGdal command is not properly processing the HDF file into the processed images, or why the gdalinfo command is not found.


texts:
I read about WebGL Earth that it uses under the hood the great CesiumJS project for the rendering of the data.

But exactly, for a developer, what are main differences between CesiumJS and WebGL Earth ?


texts: Shapefile straight line feature appearing as arcs at different zoom levels

texts:
I have a Shapefile that contains sanitary sewer pipe sections. ArcGIS and AutoCAD Map both display all of the pipe segments as straight lines. However, when displayed in QGIS, depending on the zoom level, some of the pipe segments will appear as arcs.

Pipe as Arc

Pipe as Straight Line

What do I have to do to fix this?


texts: How to get radiobutton checked after pop-out window?

texts:
I am trying to migrate my program from ArcMap 10 VBA to ArcObjects VB.net add-in in Visual studio express 2008.
I write a ArcGIS addin component tool to draw a geometry, perform a spatial query and then pop out a windows form with radio buttons to decide based on which layers in TOC to perform the query. Now what I have encountered is in OnMouseDown Event the radiobutton read from current form and return no value. when I execute the code it opens the form as i expected but doesn't allow me to select radiobutton, besides it doesnot show any error or relevant warnings in VS express errorlist.

Here are the codes in Tool Class

Public Class Tool1
Inherits ESRI.ArcGIS.Desktop.AddIns.Tool
Dim frm2 As New Form2

Public Sub New()

End Sub

Protected Overrides Sub OnUpdate()
    Enabled = My.ArcMap.Application IsNot Nothing

End Sub

Protected Overrides Sub OnMouseDown(ByVal arg As ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs)

    ' MyBase.OnMouseDown(arg)
    Dim pMxDoc As IMxDocument
    Dim pMap As IMap
    Dim pActiveView As IActiveView
    pMxDoc = My.ArcMap.Document
    pActiveView = pMxDoc.FocusMap
    pMap = pMxDoc.FocusMap

        Call AddDataFromFile.AddShapeFile("\\Hkhkgsv02\GIS\Basemap\Database_Index", "B1k_index")
        Call AddDataFromFile.AddShapeFile("\\Hkhkgsv02\GIS\Basemap\Database_Index", "B5k_index")
        Call AddDataFromFile.AddShapeFile("\\Hkhkgsv02\GIS\Basemap\Database_Index", "B10k_index")
        Call AddDataFromFile.AddShapeFile("\\Hkhkgsv02\GIS\Basemap\Database_Index", "B20k_index")

    'Determine Selection Layer based on which scale
    Dim pScreenDisplay As IScreenDisplay
    Dim pRubberEnv As IRubberBand

    pScreenDisplay = pActiveView.ScreenDisplay
    pRubberEnv = New RubberEnvelope

    Dim pEnvelope As IEnvelope
    pEnvelope = pRubberEnv.TrackNew(pScreenDisplay, Nothing)
    'Rectangle Display

    If pEnvelope.IsEmpty = True Then Exit Sub
    With pScreenDisplay
        .StartDrawing(pScreenDisplay.hDC, ESRI.ArcGIS.Display.esriScreenCache.esriNoScreenCache)
        .FinishDrawing()
    End With
    Dim pGeometry As IGeometry
    pGeometry = pEnvelope

    'Build a Spatial Query Filter

    Dim pSFilter As ISpatialFilter
    pSFilter = New SpatialFilter
    pSFilter.Geometry = pGeometry
    pSFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelContains
        frm2.Show()
         If frm2.RadioButton1.Checked = True Then
           pFlayer = pMxDoc.FocusMap.Layer(FindLayerIndex.Find_Layer_Index("B1k"))
        ElseIf frm2.RadioButton2.Checked = True Then
            pFlayer = pMxDoc.FocusMap.Layer(FindLayerIndex.Find_Layer_Index("B5k"))
        ElseIf frm2.RadioButton3.Checked = True Then
            pFlayer = pMxDoc.FocusMap.Layer(FindLayerIndex.Find_Layer_Index("B10k"))
        ElseIf frm2.RadioButton3.Checked = True Then
            pFlayer = pMxDoc.FocusMap.Layer(FindLayerIndex.Find_Layer_Index("B20k"))
        End If
     Dim pSelectionLayer As IFeatureSelection
    pSelectionLayer = pFlayer
    If pFlayer Is Nothing Then Exit Sub
    Dim pFeatureClass As IFeatureClass
    Dim pFCursor As IFeatureCursor
    pFeatureClass = pFlayer.FeatureClass
    pFCursor = pFeatureClass.Search(pSFilter, True)

    ' Set Array to Store Feature Value
    Dim NameValueArray() As String
    ReDim NameValueArray(pFeatureClass.FeatureCount(pSFilter))

    'Find Fieldindex to read Tile's Map Name

    Dim intFieldIndex As Integer
    intFieldIndex = pFeatureClass.FindField("TILE_NAME")
...

End Sub

End Class

Hope I can get some clue.thanks.

Jing


texts: Using a MODIS code for DAYMET and I don't get any output

texts:
I am using the code below that was originally provided by @Kuik, and I have since converted the bands to DAYMET 'tmax' and 'tmin' bands.

Problem is that when I run the code, I receive no output. What am I doing wrong?

var dataset = ee.ImageCollection('NASA/ORNL/DAYMET_V3')
                  .filter(ee.Filter.date('2015-04-01', '2015-05-20'));
var landSurfaceTemperature = dataset.select(['tmin', 'tmax']);

// Create a geometry representing an export region.
var features = ee.Geometry.Rectangle([-85.1417893413635, 45.31413490213395, -86.125065708551, 44.65070625463291]);

// ######### ADDED CODE ############

// According to: https://gis.stackexchange.com/questions/307548/getting-temperature-data-of-given-point-using-modis-lst-data
// map over the image collection and use server side functions
var tempToFahr = landSurfaceTemperature.map(function(image){
  var props = image.toDictionary(image.propertyNames());
  var Fahr = (image.multiply(0.02).subtract(273.15)).multiply(1.8).add(32);
  // Mask where one of Night or day temperature has no data
  var FahrMasks = Fahr.updateMask(Fahr.select('tmax')).updateMask(Fahr.select('tmin'));
  // we assume that the night temperature is the min temp, and the day temperature is the max temperature
  var meanFahr = FahrMasks.reduce('mean').rename('meanTemp');
  // Calculate the GGD and make all the negative values 0 (see: https://en.wikipedia.org/wiki/Growing_degree-day)
  var GDD = meanFahr.subtract(50).rename('GDD');
  var GDDnonNeg = GDD.where(GDD.lt(0), 0).rename('GDDnonNeg');
  return ee.Image(Fahr.addBands([meanFahr, GDD, GDDnonNeg]).setMulti(props));
});

// calculate the sum of GGD values
var summed = tempToFahr.select('GDDnonNeg').sum().rename('summedGDD');

// ######### END ADDED CODE ############

var landSurfaceTemperatureVis = {
  min: 0,
  max: 100,
  bands: ['LST_Day_1km'],
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};
Map.setCenter(-85.60371794450282,44.73590436363271, 10);
Map.addLayer(summed, {min: 0, max: 50, palette: landSurfaceTemperatureVis.palette},'summed growing degree days', 1, 1);

// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: summed,
  description: 'imageToCOGeoTiffExample',
  scale: 1000,
  region: features,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});


texts: How do I fill a combobox in QGIS with only points vector layers?

texts:
I'm new to writing plugins in python for QGIS. I've been searching on internet for a while and wrote the code below. Even though it fills my combobox, all layers (rasters, vectors of all geometry types) get added. What am I doing wrong here?

#Populate the combo box with vector layers
    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer and layer.geometryType() == QGis.Point:
            self.dlg.mMapLayerComboBox.addItem(layer.name(), layer)


texts: How to programmatically access an Add-ins image?

texts:
From this thread:

Can you programmatically change the button image for a ESRI.ArcGIS.Desktop.AddIns.Button?

I've learned to change a command image programmatically, but... how can I programmatically get an image from the "Images" folder of the Addins, avoiding to write the physical file path?

Is there a way to quickly get access to the Addins resources?


texts: Can a USGS 6 inch orthoimage be given a higher resolution, like 1 foot, to improve GeoServer processing time?

texts: The USGS orthoimage is a 6 inch resolution quadrangle that is loaded into GeoServer as a WorldImage. GeoServer has an extremely hard time processing this compared with a 1 foot resolution quadrangle. How can this orthoimage be edited?


texts: PostGIS query “Error: Operation on mixed SRID geometries” only for certain tables

texts:
I have imported various shape files, all obtained from the geofabriks website, into QGIS using the 'PostGIS Shapefile Import/Export Manager'. A couple of them are throwing the error “Error: Operation on mixed SRID geometries”, on the below query.

I set the SRID during import, and have tried updating it using SET_SRID, st_transform, etc, but still get the error. Using ST_SRID to check the SRID returns 4326 (as expected).

Does anyone have any suggestions (I think I've tried everything on all the other questions about this)?

I also noticed that when adding these tables as a layer to QGIS, when zooming in a lot, the line covers the whole map (guessing this is related to the issue)

Files are from http://download.geofabrik.de/europe/great-britain/england.html (I'm having the issue only with Derbyshire, Nottinghamshire + Leicestershire, all others work as expected).

Could there be an issue with the shapefiles themselves?

Query - SELECT ST_CONTAINS(ST_ESTIMATEDEXTENT('tableName', 'geom'), ST_SetSRID(ST_MakePoint(lng,lat),4326)) as contains

Questions I've been looking at are:

PostGIS Geometry Query returns "Error: Operation on mixed SRID geometries" only for certain values
Operation on mixed SRID geometries

I've tried various methods of updating the SRID, all seemed to work, but the error persisted. I've tried setting the Point with ::geography::geometry, to force it to use SRID 4326.


texts:
I have a custom function that I'm trying to write (using Arcpy in ArcGIS 10.7.1), but I'm having a little trouble writing a successful function because I'm new to custom function writing.

I want the function to check out the spatial extension, calculate the NDVI, reclassify the image based on specific parameters, set the extent to environment settings, and save the output. Here is the custom function that I tried to write:

#Defining function
def ndviCALC():
    if arcpy.CheckExtension("Spatial") == "Available":
        arcpy.CheckOutExtension("Spatial")

    #Create red and nir float objects and perform map algebra
    var1 = Float(nir)
    var2 = Float(red)
    num = Minus(var1, var2)
    denom = Plus(var1, var2)
    ndvi = Divide(num, denom)

    #Reclassify ndvi output
    outRc = Reclassify(ndvi, "VALUE", ndviRemap, "NODATA")

    #Apply environments when saving final output
    resultRc = ApplyEnvironment(outRc)

    #Save resulting output
    resultRc.save(output2)

The file is saved as "calcNdvi" as a ".py" file.

Also, I have previously used the sys.path.append procedure to update the path directory for importing python scripts.

Here is the script that I have written, trying to use the function:

# import required modules
import arcpy
import calcNdvi

#Import classes
from arcpy import env
from arcpy.sa import *

# set environments
env.workspace = "C:\\Users\\DemoFunctions\\tooldata"
env.overwriteOutput = True
env.extent = "Boundary.shp"
env.mask = "Boundary.shp"

# define inputs, outputs, and remap table
input = "LANDSAT8_20150609.tif"
nir = input + "\\Band_5"
red = input + "\\Band_4"
ndviRemap = arcpy.sa.RemapRange([[-1.0, 0.0, 1], [0.0, 0.2, 2], [0.2, 0.4, 3], [0.4, 0.6, 4], [0.6, 0.8, 5], [0.8, 1.0, 6]])
output1 = "LANDSAT8_20150609_NDVI.tif"
output2 = "LANDSAT8_20150609_NDVI_RC.tif"

#Use custom function
try:
    #Trying to call the custom function here
    ndviFunction = ndviCALC
    print "NDVI successfully saved to workspace directory"

except arcpy.ExecuteError:
    print arcpy.GetMessages(2)

finally:
    arcpy.CheckInExtension("Spatial")

However, I get this error:

Traceback (most recent call last):
  File "<pyshell#26>", line 1, in <module>
    ndviFunction = ndviCALC
NameError: name 'ndviCALC' is not defined

How can I write this function correctly so that my function works?


texts: Python Tools: What parameter table type is best?

texts:
I'm trying to decide on the "best" choice for working with a set of input parameters. The data are table like. (several rows and columns)

As per: http://desktop.arcgis.com/en/arcmap/latest/analyze/creating-tools/defining-parameter-data-types-in-a-python-toolbox.htm

There are several Parameter Data Types related to tables

Record Set/GPRecordSet
Value Table/GPValueTable
Table/DETable
Table View/GPTableView

There is information about Value Tables, I gather they are best suited for 2 dimensions.  But I can't find a lot of information about the other generic Table types, and why one would be chosen over the other.  They all seem to behave the same when creating a tool in a python toolbox.


texts: Dissolve connected lines until line forks

texts: I have a lot of small lines which i need to connect, but only if connection point has 2 connected lines. When there is a fork with 3 or more connected lines, i need them to stay separated. Is there a way to do this using existing tools or do I need a custom script?


texts: Find census block for street intersection with tiger geocoder

texts:
I'm geocoding lots of addresses and map them to census block. Some of addresses are street intersections, like this

select ST_X((geo).geomout)::numeric(8,6), ST_Y((geo).geomout)::numeric(8,6)
    from geocode_intersection(' 7TH ST NW ','G ST NW','DC','WASHINGTON','20001',1) as geo

  -77.021914;38.898446;

I was using the gemout field from geocode_intersection output as the input for finding census block:

select a.addid, b.tabblock_id
        from intersection_table as a, tabblock as b
        where ST_Contains(b.the_geom, a.geomout)

This only return 5 records in 19 addresses.

Take the example of first address, this failed when I used 6 digits after dot

SELECT tabblock_id
    FROM tabblock
    WHERE ST_Contains(the_geom, ST_SetSRID(ST_Point(-77.021914, 38.898446), 4269))

But it will success with 5 digits

SELECT tabblock_id
    FROM tabblock
    WHERE ST_Contains(the_geom, ST_SetSRID(ST_Point(-77.02191, 38.89845), 4269))

  '110010058002005'

Using lon and lat with 5 digits after dot will find result for all:

select a.addid, a.geomout, get_tract(a.geomout),b.tabblock_id
    from intersection_table as a, tabblock as b
    where ST_Contains(b.the_geom, ST_SetSRID(ST_Point(a.lon, a.lat), 4269))

I can understand it's difficult to have a definitive answer for street intersections because they could be the boundary, but why using 5 digits will "solve" this problem? Should I trust the result of this method?

EDIT I found something interesting:

Result from my code:

row_seq        | 2716
street_1       |  FLORIDA AVE NW
street_2       | MASSACHUSETTS AVE NW
state          | DC
city           | WASHINGTON
zip            | 20008
addid          | 21
rating         | 3
lon            | -77.04879
lat            | 38.91150
output_address | 2198 Florida Ave NW, Washington, DC 20008
geomout        | 0101000020AD100000CB9717601F4353C0BAA29410AC744340

Use the lon/lat value to look up census block by tiger geocoder:

SELECT tabblock_id
    FROM tabblock
    WHERE ST_Contains(the_geom, ST_SetSRID(ST_Point(-77.04879, 38.91150), 4269));
> '110010041003022'

Use this lon/lat value in census geocoder website give same 110010041003022.

However, if I used the street address in Census geocoder, it give the coordinates and census block as

X: -77.048775 Y: 38.91151
GEOID: 110010055001010

I guest it's because geocode_intersection got the two street intersection by calculation, then use the coordinates to reverse geocode an address, which is the closest one to the location.

However if you input an address to look up census block, it will get the address location first, which is not exactly the same with the street intersection coordinates.


texts: ArcGIS Publish query layer as Feature Service fails

texts:
I have a query layer that's simply a DB view call in PostgreSQL, the view is as such:

   CREATE OR REPLACE VIEW daily_query AS
 SELECT b.id, b.first_name, b.last_name,
    c.assignment_date::date AS assignment_date, d.case_id, d.case_sequence,
    d.status, a.address_complete, a.the_geom_point
   FROM universe a, staff b, enumerator_assignments c,
    enumerator_case_assignments d
  WHERE b.id = c.enumerator_id AND c.objectid = d.enumerator_assignments_id AND a.id = d.case_id;

The query layer is simply a call to the view

select * from daily_query

The layer renders fine in ArcMap - but then I try to publish the layer as a Service to ArcGIS Server,  I want this to be a FeatureService so I enable "Feature Access" Capability and I only allow "Query" operation.  and I don't allow "Geometry Updates".

When I run "Analyze" it returns with "00134: Layer's data source is not supported"

The Service Analyzes and Publishes fine if I don't enable "Feature Access" on it.

I'm running 10.2.2

Any ideas?


texts:
Some times in my tools I need to output Cyrillic text via arcpy.AddMessage function.

Previously in 10.0 it was very easy. I had .py file of the tool saved in UTF8 codepage.
This code wrote messages in Cyrillic with system locale 'cp1251':

source_message = u'Some Cyrillic text here'
message = source_message.encode('cp1251')
arcpy.AddMessage(message)

But in 10.1 the same piece of code writes <msg>0; to messages.

I have tried several ways of output (utf8, cp1251, unicode) with combinations of saving source .py file in different codepages - no way.

May be anybody faced similar problems in 10.1 ?

PS I am trying to run this tool in ArcCatalog 10.1 sp1


texts: How to convert images from the South African co-ordinate system to WGS84 lat/long?

texts: I have a number of orthophotos for an area in South Africa, but they are projected in the South African co-ordinate system which measures a location in meters from a central (uneven) meridian.  These particular images should use Lo31.  However, QGIS does not seem to recognise the South African co-ordinate system and I am battling to display the images correctly.  Any ideas how I can solve this dilemma?


texts: Inserting Point into PostGIS geometry column from golang?

texts:
I have this go code to insert data in to postgresql:

sqlStatement := `
    INSERT INTO my_table (
        column_name, geog)
    VALUES ($1, ST_GeomFromText('Point($2 $3)'))
`

and then

_, err = db.Exec(sqlStatement,
        'whatever',
        longitude,
        latitude)

where latitude and longitude are float64.

I get the following error:

pq: got 3 parameters but the statement requires 1

I also tried the following, using the go-postgis library:

sqlStatement := `
        INSERT INTO my_table (
            column_name, geog)
        VALUES ($1, ST_GeomFromText($2))
    `
_, err = db.Exec(sqlStatement,
            'whatever',
            postgis.PointS{4326,longitude, latitude}
)

In this case, the error message is

pq: invalid byte sequence for encoding "UTF8": 0x00

Any hints on what is the correct way to do this simple operation?


texts:
I'm trying out QField 1.4.4 and get this error when trying to attach photos to point data.
Value "DCIM/JPEG_20200303_130156.jpg" NOT NULL could not be converted to a compatible value for field Photo_1(String).
Using a Google Pixel with Android version 10.


texts: How to use the "Grid for Indexing" in the print composer QGIS 2.14.2?

texts:
I created a grid with plugin Create Grid for Indexing. Until now all fine.

But how can I use the grid (Create Grid for Indexing) in the print composer?

The task:

I have a map and a overview. I want to have information about the location of the map, for example, "The map is in M1" or "The map is in G5", etc.

How to make that in the print composer and using Create Grid for Indexing?


texts: Save as rendered image does NOT save in CRS/Projection resolution

texts:
I am trying to save a QGIS layer (the layer in question is an STRM geotiff, downloaded from Earthexplorer with resolution 1801 x 3601) in the project CRS projection (WGS 84 EPSG:3857), as an image file. However, no matter which options I choose it ALWAYS saves in the native resolution of 1801 x 3601. The layer appears perfect in the QGIS project window, but the saved-as image is always 1801 x 3601. (the WGS 84:3857 should be approx 1801 x 3200 I think, based off a basic measurement of the width x height in the project window). In the save as options, there is even a drop down list to select the CRS, but it has no effect on the final image produced.

(The screen grab below is not using 3857 CRS but still demonstrates the problem I am having)

I can change the resolution manually, I know, but for my project I need the rendered image to be exactly as per the WGS84 EPSG:3587 resolution.


texts: Using ST_Transform when inserting new row?

texts:
I have written a script which reads a CSV file with lat/lon coordinates in WGS84 and inserts them into a table in PostGIS.

What I want to do now is to change the insert query in order to make a transformation of the coordinates from WGS84 to another reference system.
My qyestion is if something like this is possible while inserting the data.

I was trying something like this but doesn't work:

   INSERT INTO imgRural(geom) VALUES (ST_Transform(ST_GeomFromText('POINT(23.771111111111 38.013888888889)', 2100)),2100) RETURNING imgId;

What am I missing here?


texts: Visualize geographic spatial data from SQL Server directly in PowerBI

texts:
Is it possible to directly query geographic spatial data form SQL Server in PowerBI like its possible with Tableau?
https://cloudblogs.microsoft.com/sqlserver/2018/05/14/visualize-geographic-spatial-data-from-sql-server-directly-in-tableau-2018-1/

I know there is the "workaround" with ArcGIS Online but this means I have to keep two datasets up to date, the ArcGIS Online shapefile and the spatial data in SQL Server.


texts: Depicting Table of Imported (CSV) Points With Lat/Long in Another Table in QGIS

texts:
I have a CSV table of points with descriptive location. Say something like this:
Sarah, Todtmoos
Henry, Paris
John, New York
Elma, Paris
... more rows like this ...

I have another CSV table with the geographic location of these (often repeated) locations. Say something like:
Todtmoos,47.739348, 8.001628
Paris,48.853868, 2.343282
New York,40.737401, -73.994330

I would like to depict the points in the former table as located in the appropriate location which can be found in the latter table (obviously points will overlap due to duplicate locations).
I could use LOOKUP in a spreadsheet application like Excel before I bring the material into QGIS through CSV import with geography from lat/long, or, more tediously, simply do a find-and-replace on each location and replace with its equivalent lat/long.
However, I've been wondering if there is a simple way to do this within QGIS. Something with a field calculator or some other dynamic lookup process with the attribute tables? I'm guessing this is not a normal table or spatial join and would welcome suggestions on most robust way to do this. (I use QGIS 3.10 on OSX). I have seen something about conversion to a database but am hoping for something simpler?


texts: Google maps API charset for http request ?? Geocoding Requests

texts:
im sending the following request to maps api but receiving the json response with chars problem. Can some one point me in the right direction to set charset correctly ?
Request -- http://maps.googleapis.com/maps/api/geocode/json?latlng=38.676817,-9.178433&sensor=false
Response of formated_address = R. dos LusÃadas, 2805 Almada, Portugal
Where it should be R dos Lusíadas, 2805 Almada, Portugal.

Thanks
Jorge


texts: Google Cloud Image exporting blurry and small tiff images

texts:
I'm trying to export a image from GEE. For some reason, I'm getting a very small image size that is only 21x30 in dimensions.
I'm trying to read it in Python and QGIS. Same problem in both platforms
GEE Code:
Export.image.toDrive({
  image: filtered_region,
  folder:"pixxel",
  region: geometry,
  scale:30,
  fileFormat: "GeoTIFF"
})

Link: https://code.earthengine.google.com/b79441422defdf9a7ec08ca75cedccbf
Python Visual Output:

Why is it getting exported in this way?


texts: Adding label to map layer by clicking feature or similar action in QGIS?

texts:
I have searched the web (including stack exchange) quite extensively and haven't found anything. I've tweaked the 'repeat' setting on the label position tab on a line layer to where it labels enough features without me having to delete a ton of duplicates. But there are still some lines that are not labeled at all (complex stream network). I'd like to add a label for these streams. Ideally, I could click the stream and it would populate the name into the label automatically. But I'd be plenty happy to click the label location and type it in by hand. I just need to get the label on the map.

Am I missing something?

I am using QGIS 3.x.


texts: Ways to measure similarities/congruence of tracks (point features) in ArcGIS

texts:
I have a sets of GPS coordinates recorded while tracking different cross country runners. What I would like to do is to compare how similar the tracks of the individual runners are.

Is there a test for congruence or some such in ArcGIS 10?


texts: Create Raster file in R from Matrix and latlons

texts:
I could use some advice on how to create a Raster file in R with a Matrix of data points and the coordinates stored in a separate file.
Both, data and coordinates come from a NetCDF file. I extracted the data and made some calculations that resulted in a matrix looking similar to this one:

NA NA NA 12 NA NA NA
NA 23 34 13 45 NA NA
23 12 12 23 34 54 56
NA NA NA NA 23 21 NA

With the difference, that there are 79 rows and 78 columns.
The unprojected coordinates represent the centroids of the raster and can be extracted from the NetCDF file. Logically the longitude coordinates have a length of 79 and the latitude of 78. i saw a very good answer on a similar problem here on GIS Stack Exchange but I'm not quite sure if I can apply the solution to my problem since my data structure is different and my coordinates are not projected. Maybe one solution would be to project the coordinates and then make a long table reshaping the data with the coordinates. But maybe there is also an easier solution...

csv files containing the data and the original netCDF can be found here. Any help would be very appreciated!!!


texts: Calculate field interactively in ArcGIS ModelBuilder

texts:
I'm working on a closest facility analysis using ArcGIS ModelBuilder and want to enable the user to add locations interactively (i.e. typing the coordinates, not using any table or source layer).

I created a point feature class with two fields (X, Y) in ModelBuilder and added a new row -a new point- using a small Python script, so I only need to create something that allows the user to type the coordinates of the point and copy them to the (X, Y) fields of the feature class.

I have tried it with Create Variable (either Point or Double) + Collect Values + Append, but it does not work because the outputs of Collect Values are not recognised as valid inputs in Append.

Any suggestions?


texts: Data Driven Pages- Can i grey out non-zoomed to layers of each page?

texts: I am making a map book of plots in a municipality. I am using data driven pages, but I want each page to only show the current plat and grey out the other plats on the edges. Is there a code I can use to do this? Using 10.1


texts:
I wish to convert files (1000s) from NetCDF4 format to geotiff using an arcpy script. When I try to run the script, I do not see any output and after the fourth file, the python script abborts. Is this a problem with my code? Here is the code:

import arcpy
from arcpy import env
import os,sys
from arcpy.sa import *

arcpy.env.overwriteOutput = 1
arcpy.CheckOutExtension("Spatial")

inPath=r'E:\TRRM\daily'
outPath=r'E:\TRRM\output'
env.workspace = inPath
arcpy.env.scratchWorkspace = inPath
nc4List = arcpy.ListFiles('*.NC4')

try:

    var = "precipitation"
    xDim = "lon"
    yDim = "lat"

    for nc4 in nc4List:
        inNetCDFFile = nc4
        outRasterLayer = os.path.join(outPath, nc4.split(".")[1] + ".tif")
        arcpy.MakeNetCDFRasterLayer_md(inNetCDFFile, var, xDim, yDim,outRasterLayer,"","", "BY_VALUE")
        print "Making RasterLayer from ....."+str(nc4)
    print "Done"

except:
    print arcpy.GetMessages()


texts: Insert a MultiPolygon in a PostGIS table with MapInfo?

texts:
I'm using MapInfo 16 to edit a PostGIS table through ODBC driver.
The geometry column type of my table is MultiPolygon.
The value in spatialtype colmun of a record in mapinfo.mapinfo_mapcatalog table is 19.2. That stands for MapInfo region type AFAIK.

When I'm hitting a Save Table button, I get an ODBC Error:

Geometry type (Polygon) does not match column type (MultiPolygon);

So I guess MapInfo is trying to insert a Polygon geometry into a MultiPolygon column. Converting new created features to regions doesn't help.

The only workaround I can think of is to set a BEFORE INSERT trigger on this table with something like

NEW.geom := st_multi(NEW.geom);

But it doesn't seem like a natural solution.
Maybe I'm missing something in MapInfo options?

UPDATE:Proposed trigger solution doesn't work. Seems like type checking happens earlier than trigger fires.


texts: Adjust polygon fill transparency via ArcGIS cartographic representations?

texts:
This is my first experience with ArcGIS representations (10.3.1, advanced).  I would like to adjust polygon fill transparency by editing the representation (via Layer Properties>Symbology>Representations).  Although I can modify the color pattern, symbol levels, etc. via representations, I don't see the option to change fill transparency. Here's a screen shot:

I've reviewed a 2010 ESRI pdf (cartographic-representations-tutorial.pdf), but transparency is not mentioned, so I've resorted to using the Layer Properties>Display>Transparent setting.  Shouldn't I be able to adjust transparency via representations?


texts: R kriging raster returns NA for some prediction points?

texts:
I am trying to interpolate a raster layer of chlorophyll-a data using kriging (gstat::krige) in R, however the resulting interpolation is showing NA values in a peculiar straight line through the raster (see last plot).

From the interpolated raster I have then extracted values corresponding to spatial points.

  The chlorophyll-a raster (chl), colour palette (pal) and spatial points (subSP) relevant to this question can be found here.

library(raster)
library(gstat)
library(rworldmap)
data(countriesLow)

#Load subsetted lat and lon coordinates
load("track_points.Rdata")
head(subSP)
coordinates(subSP) <- c("lon", "lat")
projection(subSP) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

#set grid xmin, xmax, ymin, ymax of total dataset extent of ALL track points
parGrid <- raster(extent(-72, 52, -68, -21), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
res(parGrid) <-  1.851111 ## ~about 100km?? Not sure if this is correct.
#Crop the parent grid to the extent of the SUBSET of track points (plus a little buffer)
ex <- extent(plyr::round_any(extent(subSP)[1]-2, 1, floor), plyr::round_any(extent(subSP)[2]+2, 1, ceiling), plyr::round_any(extent(subSP)[3]-2, 1, floor), plyr::round_any(extent(subSP)[4]+2, 1, ceiling))
subGrid <- crop(parGrid, ex)

#Read in chlorophyll-a raster (same extent as cropped grid).
load("chl_raster.Rdata")

#Plot raster to identify missing values
load("chl_pal.Rdata")
plot(chl, col=pal$cols, breaks = pal$breaks, legend = FALSE, main="Chl-a", xlab="lat", ylab="lon")
mp <- crop(countriesLow, extent(chl))
plot(mp,col="black",border=FALSE, add=T)
plot(subSP,add=T, pch=16, cex=0.3, col="grey50")

#Krige raster layer to interpolate NA values
#Convert chl raster to dataframe
P <- data.frame(X=coordinates(chl)[,1], Y=coordinates(chl)[,2], chl=values(chl))
#Remove NAs from chl raster
P <- P[!is.na(P$chl),]
#Convert to SpatialPoints
coordinates(P) <- c("X", "Y")
#Set projection
projection(P) <- projection(chl)

#Convert empty raster to SpatialGrid
grd <- as.data.frame(coordinates(subGrid))
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE  # Create SpatialPixel object
fullgrid(grd) <- TRUE  # Create SpatialGrid object
proj4string(grd) <- proj4string(P)

# Replace chl point boundary extent with that of empty grid - perhaps this is where the problem orginates??
P@bbox <- grd@bbox

# Define the 1st order polynomial equation
f.1 <- as.formula(chl ~ X + Y)
# Compute the sample variogram
var.smpl <- variogram(f.1, P)
# Compute the variogram model by passing the nugget, sill and range values.
dat.fit  <- fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE, vgm(psill=0.65, model="Exp", range=800, nugget=0.15))
# Plot variogram
plot(var.smpl, dat.fit, xlim=range(var.smpl$dist), ylim=range(var.smpl$gamma))

# Perform the krige interpolation
dat.krg <- krige(f.1, P, grd, dat.fit)
# Convert kriged surface to a raster object
r <- raster(dat.krg)

#Plot interpolated data
plot(r, col=pal$cols, breaks = pal$breaks, legend = FALSE, main="Chl-a interpolated", xlab="lat", ylab="lon")
plot(mp,col="black",border=FALSE, add=T)
plot(subSP,add=T, pch=16, cex=0.3, col="grey50")

What is causing the straight line of NA values on the left side of the raster plot and how do I fix this?

#Extract chl from interpolated raster values
chl_ex <- raster::extract(r, subSP)


texts:
The purpose is to create a small web application using an open source product, preferred Leaflet (but I will Switch the product if it works within another Framework than leaflet), simply showing a base map and with the possibility that the user can load (through a search on the computer system) raster images that only have to be displayed.

What is the basic conceputal Approach to solve this problem?


texts: Detecting building footprints orthogonality using ArcMap

texts:
I have a number of buildings vector data from OSM and I am trying to automatically differentiate those buildings which aren't orthogonal in its edges from those that are and then calculate the error angle +/- that doesn't make it orthogonal. The error angle values will be saved as a new field (those with true orthogonality get an error angle of 0.0 as its field value).

Is there a way to do this in ArcMap 10?


texts: Web page that dynamically draws contours at user specified elevation?

texts:
How might I setup a web page that will draw dynamically contours and/or hypsometric tints at a user specified starting elevation and/or interval?

Example scenario: flood hazard modelling. The web page supplies the elevation model and widget controls, the user asks the page "what would be underwater at n meters?" Assume the horizontal extent of elevation model would be small, less than 2000x2000 cells/pixels.

Client side processing preferred. I know to how to set up a processing pipeline in ArcGIS Desktop and QGIS to accomplish this (derive contours and contour plugin respectively) but I want to take myself out of the middle and let people do this for themselves. The interaction pattern of "input parameters >> send for processing >> wait for results >> 'oh, what about N instead?' >> ... " has too much friction.

ArcGIS Server and ArcGIS Online are available, Mapbox is an option, but a self contained setup that could run on any reasonably modern web environment is preferred.


texts:
How can I draw Distributions Maps (double circles) in ArcGIS 10.1 ?

like this:

This map drawn in CorelDraw,it is shows the number of schools and the number of students in each Province (the circle with Vertical lines is the number of students,the circle with
Horizontal lines is the number of schools).

So,I want to make this circles in ArcGIS (like the Divided circles in ArcGIS )


texts: Cell Size parameter for IDW (Spatial Analyst) tool

texts:
I am attempting to use the IDW (Spatial Analyst) tool in a Python script however I am a bit confused on what to put for the "Cell Size" optional parameter. When using the tool within ArcMap it seems to calculate this variable for you when the input is selected(pictured),however I am unsure if when scripting I would have to calculate it manually or if I leave it blank ArcMap will calculate it when I run the script?

Not sure if it matters but the script is not Stand Alone it will be used as a scripting tool run in ArcMap


texts: Where is the licence manager installed on ArcGIS Server 10.3?

texts:
According to the ESRI Documentation, the licence manager for 10.3 aka. LSAdmin.exe is located here:

<installation_path>/arcgis/license10.3/LSAdmin

On my machine it is not there. I want to install and licence the image server extention but this is not going to happen if I can´t authorize it.


texts: I'm trying to add high resolution orthoimagery from the USGS to a mosaic dataset in arcMap. When I "Add rasters to mosaic dataset" most of the rasters load in fine, but some have strips of white area and other places are sort of grey, low contrast. I checked Update Overviews. Please see the image below as an example of the issue. Also, as you zoom in the grey/white areas go away. Any ideas?