texts: How to select polygons associated with (not necessarily containing) corresponding points?

texts: I have a list of polygons, each associated with the x/y coordinates of a point. These points may not intersect the polygons. How can I select the polygons that do contain the corresponding point? I only have a ArcView license.


texts: What does it mean for two datums to be navigationally equivalent?

texts:
I read from NOAA that NAD83 and WGS84 are "navigationally equivalent"

  NOAA uses the North American Datum of 1983 (NAD 83) for all of its
  charts. This datum is navigationally equivalent to the World Geodetic
  System of 1984 (WGS 84).

I am not familiar with this concept. Could some one explain how the two datums are different, and how are they equivalent in terms of navigation?


texts: Checking if layers in dataframe using ArcPy?

texts:
I working on mxd (I work with ArcMap 10.2.2 & python 2.7.5.) that have 100 layers. I'm trying with python script to add layers to dataframe and check if the layers are within the dataframe. if it within I want that python will add the layers to the table of content, and if it not within the dataframe- the python will remove the layers. Manually I can do it by adding each layer and see if it in the dataframe and then start to choose what to remove, but it will take a lot of time.
I also can use this option:

and it also will take a lot of time. I have this script:

import arcpy
from arcpy import env
area= '162000 631000 172000 641000'
env.workspace = r"C:\project"
for mxdname in arcpy.ListFiles("*.mxd"):
print mxdname
mxd = arcpy.mapping.MapDocument(r"C:\project\\" + mxdname)
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
df.extent = area
for lyr in arcpy.mapping.ListLayers(mxd, "",df):
    if lyr.name == "residence":
        addLayer = arcpy.mapping.Layer(r"C:\project\layers\residence.lyr")
        arcpy.mapping.AddLayerToGroup(df, lyr, addLayer, "BOTTOM")
        else lyr.dataSource == r"F:\GIS\topo_50000\50000.sid":
            arcpy.mapping.RemoveLayer(df, lyr)

    if lyr.name == "rivers2":
        addLayer = arcpy.mapping.Layer(r"C:\project\layers\rivers2.lyr")
        arcpy.mapping.AddLayerToGroup(df, lyr, addLayer, "BOTTOM")
mxd.save()
del mxd

I hope someone can describe an easy way to accomplish what I'm after with python script?


texts: Calculate average multilinestring point elevation

texts:
I'm triying to calculate the average elevation of some multilinestring points.

Is there a fastest way than this pgsql function?

CREATE OR REPLACE FUNCTION topstation.sig_moyenne_z(ingeom geometry)
  RETURNS double precision AS
$BODY$
DECLARE
    geom2 geometry;
    sum double precision;
    max double precision;
BEGIN
    sum:=0;
    geom2=ST_LineMerge(inGeom);
    max:=ST_NPoints(geom2);
    IF max > 1 THEN

    FOR i IN 1..max
    LOOP
        sum := sum + ST_Z(ST_PointN(geom2,i));
    END LOOP;
    RETURN sum / max;
     ELSE
    return 9999;
     END IF;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE
  COST 100;


texts: Coordinate array to layer using cheap-ruler

texts:
I have calculated a destination from bearing/distance (feet) using cheap-ruler.js https://github.com/mapbox/cheap-ruler

Console.log() returns this array from the calculation..

[-87.36328091530362, 37.78808138412046]

I then am trying to add the result to the map using...

    var ruler = cheapRuler(35.2, 'feet');

    var bearing = document.getElementById('bearing').value;
    var distance = document.getElementById('distance').value;

    var calc = ruler.destination(lnglat, 0.1, 90);

    destlnglat = [calc[0], calc[1]];

    var myIcon = L.divIcon({className: 'css-icon'});

    // Pass array to layer...
    var destinationPoint = L.geoJson(destlnglat, {
    pointToLayer: function (feature, destlnglat) {
            return new L.Marker(destlnglat, {icon: myIcon });
        }
    }).addTo(map);

Which returns no error in the console? How do I correctly pass back to L.geoJson()?

Edit: I got this working using the following..

function calculateDestinationFt() {

    var ruler = cheapRuler(35.2, 'feet');

    var bearing = document.getElementById('bearing').value;
    var distance = document.getElementById('distance').value;

    var calc = ruler.destination(lnglat, bearing, distance);

    var myIcon = L.divIcon({className: 'css-icon'});

    var m_coords = {
        lng: calc[0],
        lat: calc[1]
    };

    var destination = {
                        "type": "FeatureCollection",
                        "features": [{
                            "type": "Feature",
                            "properties": {

                            },
                            "geometry": {
                                "type": "Point",
                                "coordinates": [m_coords.lng,m_coords.lat]
                            }
                        }]
                        }

    var destinationPoint = L.geoJson(destination, {
    pointToLayer: function (feature, lnglat) {
            return new L.Marker(lnglat, {icon: myIcon });
        }
    }).addTo(map);
}

But my results are waaaay off? I got these results pluging in a bearing:270 and distance 5280ft? Neither ft or bearing appear correct. AaAH!!


texts: I am trying to open a .img image fil in ERDAS Imagine and I keep getting the message "End of file reading dictionary" with ERDAS not responding. The program gets stuck here and I have to shut it down. How can I view the image?


texts: MongoDB: Search for polygons which have a certain location

texts:
Good day :).

I am trying to implement a search feature on MongoDB using Mongoose.js , the short version of which is this:

I tag locations on a map as polygons. and give them names. Then when an user searches for a place using coordinates, I check which polygon the location falls within, and return the polygon.

I know that the reverse of what I am trying to do works by specifying a polygon, and then searching for locations which fall within the polygon, but is the opposite possible?

So instead of saying:

var polyA = [ [ 10, 20 ], [ 10, 40 ], [ 30, 40 ], [ 30, 20 ] ]
query.where('loc').within.polygon(polyA)

I want to say

var 'loc' = [ 11, 24 ]
query.where(loc).within.polygon('polygon')

where 'polygon' is a polygon schema.

I believe it might be possible to do this with mongoDB 2.4 as it'll support GeoJSON, but is it currently possible to do it with mongoDB 2.2?

Appreciate the help.

What I've thought of so far

I think it might not be possible even with MongoDB 2.4. So I'm thinking of using some linear algebra to solve this.

I compute the center of the polygon, and store it with the polygon area.
Then when I search the database I find a polygon with the center/area-point within a certain distance of the location.
Loop through each two points in the polygon, and set an y <?> m(x) + c type constraint, store the constraints in a JavaScript object.
Check if location coordinates satisfy the constraints, and return false if any of them is not satisfied. Else return true if satisfied.

In effect I'm drawing a polygon's restrictions, and checking if the said location meets all of these. The problem here's that if I have complex polygons (which I'll more likely have) the test would take longer to return an answer.

On reinventing the wheel

I sometimes get reprimanded for 'reinventing the wheel' by using MongoDB instead of a proper GIS database, which normally has some of the features I ask about. I am aware that what I am trying to do is what has probably been solved years ago by RMDBS, but I want to do it in Mongo.
One prevalent way which this feature is likely already implemented is in geocoding, because to find the town/state/country which a point lies in, one would do the above. I thus sincerely want to try to do this in MongoDB as it complements the other computations that I'm already performing in Mongo.
Thanks again


texts:
I am using what I would describe as the most primitive method of doing intersection with OGR. The short script below describes how I do it. Is there a best way of doing this?

   from osgeo import ogr, osr
   shp1 = ogr.Open(file1)
   shp2 = ogr.Open(file2)

   SpatialRef = osr.SpatialReference()
   SpatialRef.SetWellKnownGeogCS('WGS84')
   # Create dst file here
   dstshp = ogr.CreateDataSource('SomeFilename.shp')
   dstlayer = dstshp.CreateLayer()
   # define its attribute fields for dstlayer and create them

   layer1 = shp1.GetLayer(0)
   layer2 = shp2.GetLayer(0)

   for feature1 in layer1:
      geom1 = feature1.GetGeometryRef()
      attribute1 = feature1.GetField('FieldName1')
      for feature2 in layer2:
         geom2 = feature2.GetGeometryRef()
         attribute2 = feature2.GetField('FieldName2')
         intersection = geom2.intersection(geom1)
         dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
         dstfeature.SetGeometry(intersection)
         dstfeature.setField(attribute1)
         dstfeature.setField(attribute2)
         dstfeature.Destroy() # and other features must be destroyed too


texts: Cartopy unable to find tilematrix for projection

texts:
I am currently trying to run the example script from cartopy
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.shapereader import FionaReader
import os
import fiona

def main():
    url = 'http://localhost:8080/wmts/1.0.0/WMTSCapabilities.xml'
    layer = 'LayerName'
    print(os.getcwd())

    fig = plt.figure()

    ax = fig.add_subplot(1, 1, 1, projection=ccrs.epsg('3857'))

    ax.add_wmts(url, layer)
    #ax.set_extent([-15, 25, 35, 60], crs=ccrs.PlateCarree())

    ax.set_title('Suomi NPP Earth at night April/October 2012')
    plt.show()

if __name__ == '__main__':
    main()

As a server I use mapproxy with wmts. The Error message is :
ValueError: Unable to find tile matrix for projection.
    Projection: _EPSGProjection(3857)
    Available tile CRS URNs:
        EPSG:3857
        EPSG:900913

I found a similar problem with no answer yet : https://github.com/SciTools/cartopy/issues/1330


texts: I'm building a static map, and want to center text in a table I've added to the layout. How can this be accomplished in ArcGIS Desktop?


texts: Writing Raster Calculator Map Algebra expression?

texts:
I have the following map algebra expression to be used in Raster Calculator.

The "Output raster" should satisfy the following conditions:

If (IsUseAWSSites == True):
  if (InputRaster exists):
    Use "%InputRaster%" pixel value
  else:
    Use "%AnotherRaster%" pixel value
else:
  Use "%AnotherRaster%" pixel value

IsUseAWSSites is a boolean model param as variable.

InputRaster and AnotherRaster are Raster Dataset.

How do I turn this into a Con expression?


texts: Using Multiple Web Adaptors ArcGIS-Server 10.2

texts:
Currently we have a brand spanking new install of 10.2 ArcGIS-Server with the web adaptor installed on a separate server in the DMZ running HTTPS with a purchase SSL cert. We cannot bind both 80 and 443 on the IIS server as it breaks the applications using the SSL for some reason.

We need to install a second on another web server for a 3rd party app that integrates GIS web mapping into its software but cannot do HTTPS, so this will be Web Adaptor #2 using only http 80. Has anyone ever done multiple web adaptors? Are there any caveats to this?


texts: Datum Transformation from WGS84 UTM Zone 35N to ED50UTM Zone 35N with a grid shift file (.gsb file)

texts:
I would like to transform WGS84 UTM Zone 35N coordinates to ED50UTM Zone 35N coordinates with a grid shift file (.gsb file I have) in python using pyproj. I tried below code but it gives failed to load datum shift file error (code is in latlong, I also dont know how convert utm coordinates). Could you please help me?

import pyproj

ed50 = pyproj.Proj(proj=latlong +nadgrids='wgs84_ed50.gsb')

wgs84 = pyproj.Proj("EPSG:4326")

a, b = pyproj.transform(wgs84, ed50, 10, 20)

print(a,b)

Thanks in advance


texts: Render a large set of features in the map and make them interactive

texts:
We have a lot of different features (points,lines,polygons) saved in the PostgreSQL with PostGIS.

And we have to display these data in a map (we use OpenLayers 2.13 at the moment), and the different features are organized as a tree, user can check the features they are interested.

There is no problem once user the features selected by user have only a little records, we can add them in the map directly.

However when user select too many features or one feature contains too many records, it is impossible to display them all once, although we only search features inside the map extent, but think about once the map is zoomed at a small zoom level.

So we tried to use geoserver to render the map at the server side, but we need the feature rendered in the map clickable, that's to say once a feature is displayed in the map, user can click it, then we will show the detail information about this feature. But we meet two problems:

1 The features are overlapped one by another.
For a given point feature layer, there maybe some points nearby each other(the distance between them are very close), so when zooom to small level, these features will overlap with each other.

2 Make them clickabe.
I know there is a WMSGetFeature in OpenLayers, but since the data are served by different wms service, so it is not a good idea.

And I have thought to request the features for the user selected feature types, and then instead of rendering them in the map, we calculate once the mouse inside the bounding of the feature, if yes, we display the feature. I am not sure if this is acceptable.

Also I have another idea, do not use wms at all, just search the data inside the viewport of the map, but we will limit the returned count of the features, and them render them with OpenLayers API.

I want to know if you guys have meet this kind of requirement, and how do you make it.


texts: How to add a hosted Feature Layer from AGOL to my local ArcGIS Pro through Python?

texts:
I`m trying to find a way to add an online layer to my local map, through Python. I know how to grab the layer info, but don't know how to use it.

onlinelayer = gis.content.search(mylayer)[0]
itemid= onlinelayer.id

So I got the item ID,  but how I apply it? Is there something similar to addDataFromPath()?


texts: How do I stop null values from being assigned extreme color values?

texts: I have a raster surface that I kriged in Surfer and imported to Grass GIS but when I view it in QGIS, the null values outside of the interpolated surface take an extreme value and wash out the resolution in the rest of my surface.  Is there an easy way to stop the null values from being assigned off the scale values?


texts: Looping a function over a sequence of parameters with iterate()

texts:
I built a function landcov(area, date) to calculate the share of land types 0 to 3 in year "date" from MODIS12Q1 within the polygons of the "area" parameter. I want to use iterate() in a loop-like fashion to do that for many "date" parameter values. The function is set as

var landcov = function(area, date) {
  // loading MODIS and filetring ICGB classification:
  var feature = ee.Image("MODIS/006/MCD12Q1/" + date).select(['LC_Type1']);
  // first feature with class==0:
  var l = ee.List.sequence(0,17);
  var feat = feature.remap(l, ee.List.repeat(0,17).insert(0,1));

  // adding subsequent layers (impossible to loop it :/):
  var feat = feat.addBands(feature.remap(l, ee.List.repeat(0,17).insert(1,1)));
  var feat = feat.addBands(feature.remap(l, ee.List.repeat(0,17).insert(2,1)));
  var feat = feat.addBands(feature.remap(l, ee.List.repeat(0,17).insert(3,1)));

  // Map land cover area per polygon (zonal stats):
  var feat = feat.reduceRegions({
    collection: area,
    reducer: ee.Reducer.mean(),
    scale: feat.projection().nominalScale(),
  });

  var feat = feat.map(function(ft){
    return ee.Feature(ft).set('date', date)
  });

  // exporting it to .csv:
  Export.table.toDrive({
    collection: feat,
    description: "land_cover_" + date,
    fileFormat: 'CSV'
  });

}

Now, suppose I want to do that function for 2001 to 2003 for the polygon called aoi. The only way I found is to run:

// load "aoi":
var aoi = ee.FeatureCollection('ft:1S4EB6319wWW2sWQDPhDvmSBIVrD3iEmCLYB7nMM').filter(ee.Filter.eq('StateName', 'Maine')).select(['CntyFips']);

landcov(aoi, '2001_01_01');
landcov(aoi, '2002_01_01');
landcov(aoi, '2003_01_01');

How would one loop over a sequence '2001_01_01' to '2003_01_01'?


texts: Extract features from satellite images and store the features as a shapefile?

texts:
I am trying to identify boundaries, blocks etc from Satellite images which are stored as GeoTIFF files. I want to extract fields, buildings etc. by feature recognition and store that as GeoJSON/shapefile files.

Are there any opensource tools/libraries which I can implement to do that?


texts: I am preparing the map of Nepal with 75 districts using the KMZ file downloaded from global admin. I opened the KMZ file (NPL_admin3.kmz) in google earth pro. I produced the map, and also worked successfully with coloring the districts. But I could not put the name of the 75 districts in the map. In other GIS tools such as DIVA or Quantum, when we 'add label' under the 'layer' menu bar, we get the name associated to the boundary. But in google earth pro, I have been facing difficulty. Can anyone please help me?


texts:
I'm looking to calculate the approximate distance in miles from a linestring.

This post recommends using ST_Transform in conjunction with the State Plane coordinate system for my area. In my case, I need to do this for many areas across the US, not just one. Is there a simpler way to calculate mileage without a need for the state-specific SRID?

For some context, I'm calculating the distance in mileage for dog walks. That means that the calculation doesn't need to 100% accurate, but being fairly close is a must. I do have the address for each walk available, which means I could use it to get the state-specific SRID, but I'd prefer to do it an automated fashion rather than having to build a lookup system on my own.

UPDATE

I'm running Postgres v9.5 and PostGIS v2.2.1.

The coordinates column is of type GEOMETRY and is in the 4326 coordinate system.

I've started with the following query:

select
    appointment_id,
    (
        ST_Length(
            ST_Transform(
                ST_Makeline(coordinates ORDER BY captured_at),
                2877
            )
        ) / 5280
    ) as distance
from path
where appointment_id = '38f6f1e3-d865-5a57-b87b-597d3f68f405'
group by appointment_id;

As mentioned above, I'm using some code from the example in the referenced post.


texts: postgis union geometries in consecutive order


texts: Manually adding a Null value to a string field in an attribute table in ArcGIS

texts: When I'm in an editing session and I manually delete a value from a numeric field the value turns to null (<Null>). Manually deleting a value in a string field creates an empty string, not a Null. Currently I select the record and field calculate Null, which is clunkier. Is there a fast way to change a string to a Null in an attribute table?


texts:
I'm trying to help a friend retrieving QGIS data from a formatted drive.

Finding the .qgs project that was saved on the desktop is no problem, but it seems that there is no actual data located in that file. When trying to load it in QGIS the following error occurs: "error occurred while parsing element on line 1 column 1 for file C:/Users/H..."

I have browsed the directories for some time now without finding finding the actual database, and I have absolutely no experience with this software or it's structure. So I'm hoping that someone can guide me a little in the dark.


texts: Set a Vector input as a variable in QGIS Tool Script

texts:
In a QGIS script, how do you set a Vector input (via your tool's GUI) as a variable which can then be used in subsequent processing in your tool's script?

I know how to set a map Layer as a variable, for example:

 pts = QgsMapLayerRegistry.instance().mapLayersByName('civic_pts')[0]

But instead of hard-wiring my file inputs like that, I want the flexibility of setting my inputs via the tool's GUI.

I've tried:

 ##Input_civic_pts=vector

...and then tried using it later in the script, for example...

 for feat in Input_civic_pts.getFeatures():

...but this generates an error.

The same script works fine when I hard-wire the civic_pts file as above.

Obviously I don't know how to convert an Input to a scripting variable.


texts: How to create cross sections from vector point layer in QGIS?

texts:
I have a set of points (X,Y,Z) across a river. Is there a way to present them as a cross section?
I'm attaching a cross section I made in Excel: The X axis are the longtitude and the Y axis, the Z elevation data.
I want to make such a cross section in QGIS.


texts: Populating a text field based on another text field using Python parsing in ArcMap field calculator

texts: I am trying to populate a text field 'Subunit' based on another text field 'ACEC_Name'using the python parser in the field calculator (ArcMap). If the ACEC_Name = "Mojave Ground Squirrel" then I want the Subunit field to read "MGS1M." Attached is my code block for field calculator, every time I run the calculator it returns  values for the entire Subunit field (even if they were previously populated).


texts: Python alternative for ArcPy's Con tool?

texts:
I'm wondering if there is a method to perform a function similar to ArcGIS's Con tool, using open source options in Python.
I have two rasters, one that shows veloticty and one that shows detph. With ArcPy  it is easy to get a new raster:
Con(((detph >= 1) | (veloticty >= 1)|(detph*veloticty>=0.5)), 3,(Con(((detph >= 0.4) | (veloticty >= 0.4)|(detph*veloticty>=0.08)), 2,1)))

I would like to do with open source options in Python.


texts: How to add line break in the ArcGIS Pro Geoprocessing window?

texts:
I'm looking to add a line break in the ArcGIS Pro geoprocessing window for a custom .pyt tool. In Arcmap, putting "\n" into an arcpy.AddMessage() statement works fine, but this does not work in ArcGIS Pro - the messages are reported in the geoprocessing window, but there is no line break. Any idea how this can be done in ArcGIS Pro?

This doesn't work:

    arcpy.AddMessage(observers)
    arcpy.AddMessage("\n")
    arcpy.AddMessage("\n")
    arcpy.AddMessage(targets)

And neither does this:

    arcpy.AddMessage(observers)
    arcpy.AddMessage("\r\n")
    arcpy.AddMessage("\r\n")
    arcpy.AddMessage(targets)


texts:
I have a layer of lines that represent stretches of electrical conductors. These sections have a point of origin and an end point. One of its attributes is the length of each of the sections.

Some of these sections have the same point of origin and the same ending point, so when I label these lines by their length, this repeated data appears for each of the lines and the data is repeated.

I need that when labeling the layer, only the label of one of the lines appears, that is, that the distance between each of the points is represented. I have tried to do it by means of a labeling based on rules, and some function of aggregation by geometry, so that only one label is shown for each space between points.


texts: Fetch continuous section between two points on a MultiLineString using shapely

texts:
Let's say I have a MultiLineString like this:

ml_string = MultiLineString([
   [[0,0], [1,1], [2,2], [3,3]],
   [[0,0], [0,1], [0,2], [0,3]]
])

Given the following two points:

a = Point([1,1])
b = Point([3,3])

Using shapely, wow can I get the first (if any) continuous section in the MultiLineString?

In this case that would be ([1,1], [2,2], [3,3])


texts: Editing of Y coordinated of polygons in QGIS

texts:
I have recently recieved an extremely complex .dfx polygon file of a mine design in Cape Lo29.The problem appears that the mine engineer has used an old local convention of using negative Y co-ordinate values, essentially plotting the design in the mid Atlantic.

I know the "fix" is fairly simple, the sign of the y-coordinate data needs to be changed to a positive. Is this possible in QGIS (version 3.4.12)?


texts: Points and bioclim layers plot correctly in R, but not in Maxent.jar

texts:
I am trying to use Maxent.jar. This program uses a CSV file of points and ascii bioclim layers (basically, a map of the Americas with climate data) and overlays them on top of one another.

I find that my points and ascii files line up correctly when I use them in R (first image), but when I input the very same files into Maxent.jar, they do not line up correctly. In Maxent.jar plots (second image), the points (which are white in the Maxent plot) are too small compared to the map and are shifted too far north.

My original data (bioclim layers and points) had a Lat/Lon projection, and the bioclim layers were rasters. I used R to project the points and the layers to the Lambert Azimuthal Equal Area projection you can see in the images. Then I used R to convert the bioclim rasters into ascii using the following code from the rgdal library:

    writeGDAL(as(MyRaster, "SpatialGridDataFrame"), MyFilePath,
        drivername = "AAIGrid")}

Obviously, there are many potential problems with geospatial data that can result in points being plotted incorrectly onto a map. What's really baffling me here is that the points plot correctly in one program and not in another.
I even read the ascii file back into R and plotted my points onto that (rather than using my original raster), and everything still works perfectly in R. I have also confirmed that the header of the ascii file (showing the the number of rows and columns and the cell dimensions) exactly matches the dimensions of my original raster.

Does anyone have suggestions for why R would plot points onto an ascii file correctly, but another program would not?

Or has anyone had this specific problem with Maxent.jar?


texts: Creating polygon from 2 txt files (Parent ID + Coordinates) using ArcPy?

texts:
I want to create polygons based on two separate text files:

the first one (as parents) lists the Polygon ID (ID_el) and three to four Points ID (IDpt1, IDpt2, IDpt3..)
the second contains the Points ID and their xyz coordinates (RW, HW,HH).

I want for each ID in the first file to set an Array (3 to 4 rows) of the coordinates in the second file

I found different scripts to do some of the parts of the process, but I really need help to figure out how to get the Points ID and coordinates for each polygon into an Array, keeping also the attributes of the parent element (ElTyp, ID_el).

Here is the bit of code I used to extract these 2 text files:

with open(arcpy.GetParameterAsText (0),'r') as f1:
    with open(r'Nodes.txt', 'a') as f2:
        f2.write("ID RW HW HH" + "\n")
        for line in f1:
            if line.startswith("ND"):
                f2.write(line.replace("ND", "").lstrip())
with open (arcpy.GetParameterAsText (0),'r') as f1:
    with open (r'Polygonen.txt', 'a') as f3:
        f3.write("ELtyp ID_el IDpt1 IDpt2 IDpt3 IDpt4 Mat" + "\n")
        for line in f1:
            coord = " ".join(line.split() [0:-1])
            mat = " ".join(line.split() [-1:])
            if line.startswith ("E3T"):
                f3.write ((coord.lstrip())+ "  " + mat+ "\n")
            elif line.startswith ("E4Q"):
                f3.write(line.lstrip())


texts: How can I load data into SQL Server spatial using SQL Server Integration Services?

texts:
How can I load data into SQL Server spatial using SQL Server Integration Services (SSIS)?

I'm aware of other tools, like FME and ogr2ogr, but my only option is using SSIS


texts: How to use interacting PostgreSQL/PostGIS triggers in QGIS to create multiple labels on polygons?

texts:
EDIT:

The trigger functions posted are well done but there are still a few issues. I've uploaded two screen recordings:

http://workupload.com/file/ygRnYtp9

http://workupload.com/file/OBY7i9Fh

I use point geometries (table 'label_point') to place multiple labels on polygon geometries (table 'soil'). After doing layer edits I've to update the affected label_point features. Therefore I've created a trigger function.

CREATE TABLE label_point (
gid serial NOT NULL,
geom geometry(point, SRID),
label_sample varchar(255),
CONSTRAINT label_point_pkey PRIMARY KEY (gid)
);

CREATE TABLE soil (
gid serial NOT NULL,
geom geometry(polygon, SRID),
label varchar(255),
CONSTRAINT soil_pkey PRIMARY KEY (gid)
);

CREATE OR REPLACE FUNCTION sample_label()
RETURNS trigger AS $body$
    BEGIN
    IF GeometryType(NEW.geom) = 'POINT' THEN
        EXECUTE 'SELECT label FROM soil WHERE ST_Within($1, soil.geom) LIMIT 1'
        USING NEW.geom
        INTO NEW.label_sample;
        RETURN NEW;
    ELSEIF GeometryType(NEW.geom) = 'POLYGON' THEN
        EXECUTE 'UPDATE label_point SET label_sample = NULL WHERE ST_Within(label_point.geom, $1)'
        USING NEW.geom;
        RETURN NEW;
    END IF;
    END;
$body$ LANGUAGE plpgsql;

CREATE TRIGGER tg_sample_label BEFORE INSERT OR UPDATE
ON label_point FOR EACH ROW
EXECUTE PROCEDURE sample_label();

CREATE TRIGGER tg_sample_label AFTER INSERT OR UPDATE
ON soil FOR EACH ROW
EXECUTE PROCEDURE sample_label();

Unfortunately this solution suffers from two issues.

1.) When deleting a soil feature or moving a soil feature (ST_Within(label_point.geom, soil.geom) = FALSE) the sample_point features are not updated to NULL.

2.) When splitting a soil feature using the QGIS 'Split Feature Tool' and changing the label of one polygon part the label_point features are not updated proberly after saving the edits.

Can anyone help me with this?


texts: Rotating symbols per attribute value in QGIS?

texts:
I am trying to apply a rule to rotate symbols based on the value of an attribute.

I tried to copy Automatically changing the angle of a symbol in QGIS?, which is essentially my question, however it does not work for me.

I have a layer with two attributes, the slope and aspect.  I want to colour the symbols by the slope (this part works) and rotate by the aspect.

The steps I take are:

Layer Properties -> Style -> Graduated -> Symbol -> Change
Select Simple Marker (arrow) and press the button to the right of angle -> field type (double) -> Aspect-sel

Applying these commands changes the symbol to an arrow, and I can easily apply the colour gradient, but the symbols do not rotate. The attribute Aspect-sel varies from 2.87 to 334.23.

ETA: I was using Version 2.8.3 - on suggestion I updated to 2.14.10, which solved the problem.


texts: How can I migrate my toolboxes and styles from ArcGIS 9.x to 10.x?

texts: I am migrating from ArcGIS 9.3.1 to ArcGIS 10 SP4. I have several toolbox and style files that I want transfer to my new computer which is now running Windows 7. Can anyone tell me where I should put them?


texts:
I have a follow up question to Statistical analyses for survey against raster?

I have a very similar problem as the initial poster David.  I am creating a multivariate model from raster data in order to predict a response variable.  In my case I am using soil and climatic data to predict crop yields.

I independently came to a very similar conclusion as the answer posted by Aaron, to use zonal stats, export to R, and create a Random Forest model.  However, I am now stuck because I want to create a raster in ArcGIS of my predicted yields based on the input rasters but I can't figure out how to translate the output from R to Arc.

In R I used this code:

library(randomForest)
fit <- randomForest(Yield ~ .,   data=data)
print(fit) # view results
importance(fit) # importance of each predictor
plot(predict(fit))

The output a plot that reasonably predicts yield.

What I want to do is take the Random Forest model and use rasters of my predictor variables (e.g. slope) from a new location, and predict the yields for the new location.  I understand how this works with simple linear regression but I can't wrap my head around how to do this with Random Forests.  What equation does RF use to make the prediction plot?  How can I get that equation and its coefficients from R to make a new map in Arc?


texts: How to identify features in an OpenStreetMap layer in QGIS 2.16?

texts:
I want identify features in a OSM layer. I saw this tool here OSM Feature widget

Where I can find it? Is still in use or is there another option?


texts:
I'm trying to join a shp and a CSV, by joinbyattibutes, by the Column "KOMNAVN"
in QGIS 2.18.

However, I get this log message error and I can't figure out what I'm doing wrong.


texts: Selecting value of zero from raster and assigning new value based on majority of its neighborhood using ArcGIS Spatial Analyst

texts:
I have a raster layer has values +1, -1, and zero. Now I want to select those pixels whose value is zero and assign them new values based on the majority of its neighbor 5*5 cells.

Is there anyway that I can do this using ArcGIS Spatial Analyst?

I've tried focal statistics but it is not exactly I want.


texts: Get list of feature types from WFS without full GetCapabilities document

texts: Is there a way to get the names of all the feature types on a WFS server, other than just calling GetCapabilities? On the server I'm interested in (GeoServer, supports WFS 1.1.0), that's a 50MB XML file.


texts: Will upgrading to Mac OS X 10.12 break my newly installed QGIS 2.12

texts: Have installed QGIS 2.12 under Mac OS x 10.9 Mavericks. Had problems when I upgraded to OS X 10.10 so reverted back to 10.9. Will QGIS stop working again if I upgrade to OS X 10.12? Anyone have any experiences?


texts: Swapping raster extents, left/right and top/bottom extents need to be reversed

texts:
I have a raster layer that was delivered with the extent boundaries (Latitude and Longitude) reversed and no spatial reference; resulting in an improperly oriented and georeferenced image:

Top: -120.5
Left: 34.0
Right: 34.2
Bottom: -120.6

What I need is:

Top: 34.2
Left: -120.6
Right: -120.5
Bottom: 34.0

My only workaround was using the rotate and flip tools to orient the image properly, but this still doesn't remedy the georeferenceing issue. I redefined the image into the proper spatial reference, but because of the extents, it still draws in the wrong location. I suspect there is a simple way to redefine the extents? Or is there a way to scrap all the erroneous parameters and input the correct ones?


texts: How to select a specific print layout (by name) in QGIS 3?

texts:
I'm a real beginner in Python. I'm using a Python script to produce maps automatically.

This script is working but it is failing in the last part: saving the current map with a specific print layout. Im using the "Maps Printer" Plugin to do this.

The code bellow worked in QGIS 2.18 but fail in Qgis 3.8.

# Your settings
composerTitle = 'Composer_001' # Name of the composer you want to export
folder = 'C:/Users/.../'
extension = '.jpeg' # Any extension supported by the plugin
titre = 'NomCarte'

mp = qgis.utils.plugins['MapsPrinter']

for composer in iface.activeComposers():
    title = composer.composerWindow().windowTitle()
    if title == composerTitle:
        mp.exportCompo( composer, folder, titre, extension )
        break

The command iface.activeComposers seems to not work anymore in QGIS 3.0.
Does anyone have an idea to solve this problem?


texts: gHydraulic plugin for QGIS 2.8 make EPANET Model Type Error

texts:
I've set up my water network layers to create an EPANET model using the gHydraulics plugin for QGIS 2.8: pipes, junctions, and a reservoir. The pipes have unique IDs, with their corresponding junctions' IDs as Node 1 and Node 2. When trying to make the EPANET model, however, I get a Python error:

TypeError: qgis._core.QgsMapLayerRegistry cannot be instantiated or sub-classed

This does not happen with QGIS 2.0.


texts:
I am trying to build a short program using cs2cs that will work on any windows pc
I copied cs2cs.exe and proj.dll to a seperate folder
When i run this to convert israel tranverse mercator to decimal latlon:

cs2cs -f "%.10f" +init=EPSG:2039 +to +init=epsg:4326

It works
but when i run the opossite:

cs2cs -If "%.10f" +init=EPSG:2039 +to +init=epsg:4326

It does not work

When i run this code from latlon in degreea to israel tranverse mercator :
       cs2cs +init=epsg:4326 +to +init=EPSG:2039

it works
what am i doing wrong?


texts: Difference between lines and polygons using geopandas

texts:
I wanted to geopandas.overlay() on lines and polygons : get the geometries from the lines that are not in the polygons.

I tried using geopandas for this :

result = gpd.overlay(my_lines, my_polygons, how='difference')

but I get this error :

  TypeError: overlay only takes GeoDataFrames with (multi)polygon
  geometries

since the overlay method only supports polygons, I did the difference using pyqgis. However I want to know if there is a way to bypass this typeError using geopandas and shapely ?


texts:
I am using QGIS 2.8.2 with Windows 7. After loading the files, I cannot edit them.
Even when I select one entity, editing buttons are still in grey and cannot be used.

Does someone know how to solve this problem?


texts:
I need help with the phrasing of this query. I am applying this in CARTO to a file with lines. The following query does not generate a curve as expected. I am new to SQL/PostGIS.

select
  mycartoid.home2015_1, the_geom
  st_transform(
    st_segmentize(
      st_makeline(origin.h_block, dests.w_block )::geography,
      10000
     )::geometry,
    3857) the_geom_webmercator
  from origin,dests

Also attached is a snippet of my data for reference:

I am trying to generate a map with curved lines similar to this link:

https://team.carto.com/u/jsanz/viz/89f6adfe-48f8-11e5-b416-0e0c41326911/public_map

SQL Query Link line 26-35 from this gist:

lines as(
  select
  dests.cartodb_id,  dests.name, dests.adm0name, dests.distance,
  st_transform(
    st_segmentize(
      st_makeline(origin.the_geom, dests.the_geom )::geography,
      10000
     )::geometry,
    3857) the_geom_webmercator
  from origin,dests


texts: Delete rows that share the same geometry value in postgis

texts:
Hello I have the following problem:

I created a grid in QGIS. For some unknown reason QGIS doesn't tesselate the plane correctly. Sometimes there are polygons with different "id" that represents the same portion of the area. I loaded the Shapefile into postgis (without noticing this) and began to have troubles in queries. When I look for the problem I saw that sometimes 2 or even 4 polygons represent exactly the same region. How could I delete those repeating polygons?


texts:
I am working on a project using Esri's Cascade Story Map app.

Within this, I am embedding web maps made using ArcGIS online. In these maps, I am trying to depict the continent of Africa, however none of the default zoom-levels in ArcGIS online are at a scale I want to use. Either they are too zoomed out, or too zoomed in. What I would like is an extent somewhere between the 3rd and 4th default zoom-levels.

Is there a way to set a custom scale to display my map at?

I have already tried changing the base map to a custom tile layer from Mapbox, however it seems that this did nothing to affect what zoom-level I can show.


texts: How to create and modify a multipolygon so that it does not overlap when merging layers in QGIS?

texts:
I am new to QGIS. I need to modify a shapefile so I can add a circle in a Canadian province to be used for Tableau afterwards.

I tried to use QGIS to modify the shapefile. I have created a new layer and I have merged it with the one containing the Canadian province boundaries. The result is as shown in the screencap below.
Instead of this result, I would like the circle to act as a province in itself so I don't have the overlapping of the province and the circle. Can this be achieved when merging 2 layers?


texts: TIN Interpolation in QGIS creates weird lines?

texts:
I am using TIN Interpolation to interpolate point features based on their cost attribute (value is between 0 and 5000). Graduated symbology creates something like this:

However when I am using TIN interpolation, I get these weird white stripes (image below). I've been trying for hours to get rid of them. They really look out of place..

Anyone have an idea how to get rid of these?


texts:
I’m having some troubles importing data from a cvs file with street addresses and the logs are not providing any meaningful information.

I installed CartoDB locally on Mac OSX 10.10.

Following are some console output which doesn't specify any errors:

Resque:

** [03:21:15 2014-11-05] 64961: Found job on imports
** [03:21:15 2014-11-05] 64961: got: (Job{imports} | Resque::ImporterJobs |  [{"job_id"=>"f248f418-5d88-4c06-ab2c-4bbbfdea88ce"}])
** [03:21:15 2014-11-05] 64961: Running before_fork hook with [(Job{imports} |  Resque::ImporterJobs | [{"job_id"=>"f248f418-5d88-4c06-ab2c-4bbbfdea88ce"}])]
** [03:21:15 2014-11-05] 64961: resque-1.23.0: Forked 65262 at 1415186475
** [03:21:15 2014-11-05] 65262: resque-1.23.0: Processing imports since 1415186475
** [03:21:15 2014-11-05] 65262: done: (Job{imports} | Resque::ImporterJobs | [{"job_id"=>"f248f418-5d88-4c06-ab2c-4bbbfdea88ce"}])
** [03:21:15 2014-11-05] 64961: Checking imports

Rails:

(0.000132s) SET DateStyle = 'ISO'
(0.000130s)  SET search_path TO "public", cartodb, public

Started GET "/api/v1/imports/f248f418-5d88-4c06-ab2c-4bbbfdea88ce" for 127.0.0.1 at 2014-11-05 03:21:15 -0800
Processing by Api::Json::ImportsController#show as JSON
 Parameters: {"id"=>"f248f418-5d88-4c06-ab2c-4bbbfdea88ce"}
 Completed 200 OK in 3ms (Views: 0.3ms | Models: 0.0ms)

 Started GET "/api/v1/imports/f248f418-5d88-4c06-ab2c-4bbbfdea88ce" for 127.0.0.1 at 2014-11-05 03:21:17 -0800
 Processing by Api::Json::ImportsController#show as JSON
 Parameters: {"id"=>"f248f418-5d88-4c06-ab2c-4bbbfdea88ce"}
 Completed 200 OK in 3ms (Views: 0.3ms | Models: 0.0ms)

The CartoDB UI shows an error:

 Unable to load data (2001)
Persisting error ? Send us the following code: "f248f418-5d88-4c06-ab2c-4bbbfdea88ce"

The csv file has following format:

TIMESTAMP,NUMBER,STREET,CITY,STATE,ZIPCODE
2014-09-05T13:16:14.178Z,800,some street,some city,some state,some zipcode


texts: arcpy.GetCount_managment not returning result in dialog window

texts:
I am having trouble with GetCount to work in my script. It is not printing the GetCount arcpy.AddMessage in the geoprocessing window dialog;

# Get the spatial reference
 spatialRef = arcpy.Describe(template).spatialReference.name

#Loop through shapfiles in folder and reproject
 for fc in fcList:
     fcspatialRef = arcpy.Describe(fc).spatialReference.name
     if fcspatialRef != spatialRef:
         arcpy.AddMessage(fc)
         arcpy.Project_management(fc, outFolder + "\\" + fc, template)
         projCount = int(arcpy.GetCount_management(fc).getOutput(0))
     else:
         arcpy.CopyFeatures_management(fc, outFolder + "\\" + fc)

# Get the count from GetCount's Result object
     arcpy.AddMessage("{0} datasets have been reprojected to: {1} and saved in {2}".format(projCount, spatialRef, OutFolder))


texts: Overlaying a matrix of intensities over the world map

texts:
I have a matrix that covers the whole world, with a different intensity per coordinate. I want to overlay these intensities on top of the world map.

I was originally using a vector layer, dividing he world coordinates by the resolution of this matrix, and creating polygons with a colour that depends on the value of this location in the matrix of intensities. When migrating it to OpenLayers 3, I have been suggested to change the implementation to using an ol.source.Raster instead.

The samples on ol.source.Raster are limited, so I am having trouble understanding the steps required. I understand I would create the ol.source.Raster more or less like:

var raster = new ol.source.Raster({
    sources: [sourceMatrix],
    operation: function(pixels, data) {
        return pixels[0];
    }
});

Where sourceMatrix has to be an ol.source.Source, so I could maybe create an ol.source.ImageCanvas with the data of the input matrix.

Would that be the right approach? Is there any other sample on ol.source.Source that would fit my use case better?


texts: Breaking a range of values into buckets for a gradient

texts:
If I have a range of values and want to make a choropleth map with a gradient, what algorithms are out there to determine the value ranges for the different colors in the gradient?

For example, if I have a dataset with the values 1,10,40,50,2000,123,23,68,39,21,10,800,600,200 and want to color 14 states based on those values, what is the best way to determine the ranges?

Something like [1-50]= light red, [51-100] red, [101-2000] dark red


texts: 'Integrating' MapXtreme with Google/Bing Maps

texts:
My company has a web application that appears to use MapXtreme in order to render it, and also uses MapMarker-USA to... retrieve the layers? We have MapInfo Pro as well, and I have no idea if that's doing anything.

Can anyone explain the differences between MapXtreme, MapMarker, and MapInfo?

Currently the application is using MapMarker 14, and MapXtreme 6.8, and I've been tasked with updating the servers running the application to use the updated versions of these programs.

One of my tasks is to improve the look of our application, and what I want to do is have the map application use either a Bing or Google map as the base layer, and our other information would appear on top of that as it has with the old base layer.

When MapXtreme 7.0 was released they were advertising 'full integration' with both Bing and Google's maps. But when looking for how to achieve this, there are no tutorials or examples that I can find. I've looked through the 7.0 release notes and developer guide, but there isn't anything relevant.

How would I go about this?

I'm really new to GIS and MapInfo products.


texts: Multiple vertices after editing giving wrong start and end x, y co-ordinates

texts:
Initially I generated some polylines across the river at various locations. I found some of them are smaller in size and not fitting the river width, Therefore, I edited those using advanced editing tool by using extend polyline tool. However,I did not have knowledge that extending would create one more vertex in the middle.

When I am computing distance using my script, the script is calculating distance between two earlier vertices it did not consider the another end vertices to compute the x,y of the actual end point of the modified lines.

How can I solve this?

I tried to delete middle vertices still I find wrong distance. I think that I may need to rebuild/repair the geometry, but I am not sure.


texts: Is there anyway to make GPModels with python scripts execute faster?

texts:
I do not know if anyone else has this problem but my GPModels that include simple python scripts seem to take way too long to execute.  One example is a model I built which is a simple select feature by attribute query and attached at the end of the model is a arcpy script that zooms to the selected feature.

Here is the code for the script.

import arcpy
MXD=arcpy.mapping.MapDocument('CURRENT')
dataframe = arcpy.mapping.ListDataFrames(MXD, "US Data")[0]
dataframe.zoomToSelectedFeatures()
arcpy.RefreshActiveView()

This model will take well over a minute from start to finish which I hope is not be the norm as this is a very simple process. Anyone have any ideas as to why this executes so slow and how to speed it up?

EDIT: The feature class it is doing the sql on is in a file GDB. In the feature class are 1100 points and each point is identified with a STORE number, there is no attribute index on. The sql expression is STORE = '(whatever store number I enter)'.  The expression seems simple enough that I doubt it could be the problem.  This GDB is on a network drive, could that be the issue, a slow network?

I just tried it again with the data local and it still took just as much time.  It must be the python script because the progress box displays "Executing Store Script" for almost the entirety of the time it takes to execute the model, the selection process takes less then a second.

I am running Windows 7 64bit if that matters at all.

I modified the script per a suggestion that the problem might be the amount of time it takes to import all of arcpy, I changed the code to whats listed below and that helped a bit.  I still think this takes too long for such a simple task.

from arcpy import mapping, RefreshActiveView
MXD=mapping.MapDocument('CURRENT')
dataframe = mapping.ListDataFrames(MXD, "US Data")[0]
dataframe.zoomToSelectedFeatures()
RefreshActiveView()

Here is what the dialog box says after the model finishes.

Executing: Model3 ""STORE" =131"
Start Time: Wed Nov 10 13:31:35 2010
Executing (Select Layer By Attribute): SelectLayerByAttribute Petm_Str_6_18_10 NEW_SELECTION ""STORE" =131"
Start Time: Wed Nov 10 13:31:35 2010
Succeeded at Wed Nov 10 13:31:35 2010 (Elapsed Time: 0.00 seconds)
Executing (ZoomToPetm): Script
Start Time: Wed Nov 10 13:31:35 2010
Running script Script...
Completed script Script...
Succeeded at Wed Nov 10 13:31:46 2010 (Elapsed Time: 11.00 seconds)
Succeeded at Wed Nov 10 13:31:46 2010 (Elapsed Time: 11.00 seconds)

Update: Exported the enitre model to a python script, cleaned it up and this is what I have now...

import arcpy

# Script arguments
Expression = arcpy.GetParameterAsText(0)
if Expression == '#' or not Expression:
    Expression = "\"STORE\" =" # provide a default value if unspecified

# Local variables:
Stores = Expression
Stores = "US Base Data\\Stores"

# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(Stores, "NEW_SELECTION", Expression)

# Process: ZoomTo
MXD=arcpy.mapping.MapDocument('CURRENT')
dataframe = arcpy.mapping.ListDataFrames(MXD, "US Data")[0]
dataframe.zoomToSelectedFeatures()
arcpy.RefreshActiveView()

Using this script only it takes 15 seconds to execute, am I being naive in thinking that this should take a couple (2) seconds to execute at most?


texts: Extracting Values from Raster image (.lyr) into Polygons using ArcGIS Pro

texts:
The raster layer is a remote sensing image (.lyr) and does not have an attribute table (cannot perform zonal stats)
Cell size is 30 by 30 meters.
I want to overlay this image over a polygon shapefile (fields/watersheds) so that I can know how many different types of pixels (area and % of crops) are located in each field.
How do I do this in ArcGIS Pro?


texts: Google Hybrid Map tiles won't load in OpenLayers when centered on certain position

texts:
Since the Google Maps API update to V3.2 (17 Feb 2015) our maps have not been loading to their initial state, centered over the Gulf of Mexico.

When the center location is: center: new OpenLayers.LonLat(-90, 25) it doesn't render, but OpenLayers.LonLat(-91, 25) works fine. Any 90 lat or 0 lon fails.

Can anyone explain why the following code does/does not work?

var map;

$(document).ready(function () {
    map = new OpenLayers.Map('olmap', {
        projection: 'EPSG:3857',
        layers: [
        new OpenLayers.Layer.Google(
            "Google Hybrid", {
            type: google.maps.MapTypeId.HYBRID,
            numZoomLevels: 20
        })],
        center: new OpenLayers.LonLat(-90, 25) //change to -91 and it works.
            .transform('EPSG:4326', 'EPSG:3857'),
        zoom: 4
    });
});

See jsFiddle for your experimenting pleasure.

If this is a bug, who do I report it to, Google or OpenLayers?


texts:
How to create geoemtry from WKT in FME?

We have a csv file in WKT format:

"ID";"GEOM"
"1";"LINESTRING ( 15, 20 16, 20 21, 12)"
"2";"LINESTRING ( 15, 20 16, 20 21, 12)"
"3";"LINESTRING ( 15, 20 16, 20 21, 12, 23, 14)"

I found  this article: https://knowledge.safe.com/articles/56486/converting-from-csv-data-2.html

But it shows how to split data into points and create a line.

Is there a simpler way?


texts: Satellites GOES16 and 17 - reprojection to Mercator is stretched

texts:
I have downloaded GEOS data from https://cdn.star.nesdis.noaa.gov/GOES16/ABI/FD/ (or https://cdn.star.nesdis.noaa.gov/GOES17/ABI/FD/). I want to reproject the disk to Mercator projection. I have used equations from https://www.cgms-info.org/documents/cgms-lrit-hrit-global-specification-(v2-8-of-30-oct-2013).pdf (section 4.4). However, resulting image is stretched.

I have set for GOES16:

"lon": -75.0,
"COFF": 5424.5,
"LOFF": 5424.5,
"CFAC": 623.331875855,
"LFAC": -623.331875855

The same data with borderd rendered. My border is red, yellow one is from GEOS data image.

Callifornia is off a bit:

However, Cuba, Haiti etc. is ok

The reprojection is correct near the center longitude and the more to the edge, the more incorrect it is.

Any ideas why? I dont have this problem with Himawari and Eumetsat using same equations (borderd are perfect fit), so I think my implementation of them is correct.


texts: Can I get gml:id from GeoServer GetFeature request

texts:
I want to make a call to GeoServer like http://someurl:8080/geoserver/wfs?service=WFS&request=GetFeature& --someParameters&...-- and get from GeoServer just items I requested in the parameters.

For example I would like to have only gml:id of all features.

Is it possible?


texts: QGIS: Interactive cursor to know latitude and longitude in boundary layer (shapefile)

texts:
Is there a way of knowing the coordinates in a specific point of a shapefile (boundary layer/map of the country) in QGIS?
Similar to what is done on Google Earth. As the cursor moves, the coordinates are changed at the bottom of the screen


texts: How to extract points off a polyline at certain distances using ArcObjects using VB.Net

texts:
I am trying to extract longitude and latitude of points from a polyline shapefile at certain intervals(10-50m). Using ArcGIS I was able to do using Densify followed by Feature Vertices to Points in Data Management Tools. But while doing so in VB.Net I am not able to locate this Class FeatureVerticesToPoints in VB.NET thus making my process useless.

Is it possible to use these ArcGIS tools in VB.NET using ArcObjects and in case yes, How?
Also, is there any other better method to do this?
Can someone help me with a code of something similar as I am a beginner to programming using ArcObjects.


texts: Getting areas of different classes in Google Earth Engine?

texts:
I would like to know if it is possible to calculate the area resulting from a classification. he been reviewing some sections, I found this

Calculating areas of different classes in Google Earth Engine?

but I can't understand the approach very well

// Traemos una imagen
var imagen = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterBounds(geometry)
.filterDate("2018-01-01","2019-01-01");

print(imagen);

var REAL = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 2000
};

var sorteo = imagen.sort('CLOUD_COVER');
var sorteada = ee.Image(sorteo.first());

var image = sorteada.clip(geometry);

// NDVI de la imagen 2018 usando expresiones
var ndsi = image.expression("(G - NIR) / (G + NIR)",{
  "G": image.select("B3"),
  "NIR": image.select("B6")
});

var pal_ndsi = {
  min: 0,
  max: 1,
  palette: ["#1F73EE","#41D3F4","#E2FA33","#FF0000"]
};

Map.addLayer(image,REAL, 'sorteada');
Map.addLayer(ndsi,pal_ndsi, 'NDSI');

Map.centerObject(geometry,10);

var feature_hielo = ee.Feature(hielo).set({class:0});
var feature_nubes = ee.Feature(nubes).set({class:1});
var feature_suelo = ee.Feature(suelo).set({class:2});
var feature_vegetacion = ee.Feature(vegetacion).set({class:3});
var feature_agua = ee.Feature(agua).set({class:4});

var fc = ee.FeatureCollection([feature_hielo,feature_nubes,feature_suelo,feature_vegetacion,feature_agua]);

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

// Con la función ee.Classifier.svm("Opciones de clasificador") se inicia el clasificador
// Es una Máquina de Soporte Vectorial (svm)
// Pero se pueden usar otras
// Las opciones están por defecto

// var classifier = ee.Classifier.svm();
// var classifier = ee.Classifier.naiveBayes();
//var classifier = ee.Classifier.minimumDistance();
var classifier = ee.Classifier.randomForest(30);

var S2Bands = ['B2','B3','B4'];

var trained = classifier.train(training, 'class', S2Bands);

var classified = image.classify(trained);

var pal = ["#FFFFFF","#AFABAB","#E2B407","#7EF03D","#6C88E6"];

var area = ee.Image.pixelArea().divide(10000);

print(area);

Map.addLayer(classified, {min: 0, max: 4,palette: pal}, 'Clasificacion');


texts: Saving QgsRasterBlock object as .tif, or to add it as QgsRasterLayer?

texts:
I am currently working on a qgis plugin and I need to simulate each pixel value in raster. I read that the best way to read each pixel value is by using QgsRasterBlock  object and I assigned simulated value to second QgsRasterBlock object object like this:

filename = "C:\Users\AB\Desktop\PyQgis\usa_dem\e10g.tif"
registry = QgsProviderRegistry.instance()
provider = registry.provider("gdal",filename)
QSettings().setValue("/Projections/defaultBehaviour", "useGlobal")
stDev = 1
raster_extent = provider.extent()
raster_width = provider.xSize()
raster_height = provider.ySize()
no_data_value = provider.srcNoDataValue(1)
block = provider.block(1,raster_extent,raster_width,raster_height)
block_2 = provider.block(1,raster_extent,raster_width,raster_height)
if block.isValid():
    for x in range(raster_width):
        for y in range(raster_height):
                        elevation = block.value(x,y)
                        if elevation == no_data_value:
                            print("no data")
                        else:
                            elevation_2 = random.uniform(elevation-
                            stDev,elevation+stDev)
                            block_2.setValue(x,y,elevation_2)

I am new to PyQGIS and I want to know is this approach is good or not.


texts:
I have a dataset of points with data at different depths:

e.g as a csv file as follows:

  "ID","DEPTH","VALUE"

  "BH01","0.5","10"

  "BH01","2.0","50"

There will be multiple points with the same ID but different depths - how can I make QGIS only label the "ID" once? Have been banging my head with this for a while and can't figure it out.


texts: Rule-based rendering of one line with two line styles

texts: I need to show parts of a network, which meet two criteria at the same time. So I would like to use the style of criteria 1 and the style of criteria 2 side by side (see picture). The problem is, that the orientation switches for every part of the network (like in the red circle). Why? I used a +0,4 offset for the green line and -0,4 for the magenta one. Where is my mistake or is it not possible?


texts: GDAL transform DTED to contour KML : gdal_contour coordinates are swapped in KML file

texts:
Using GDAL 3.1.3 (installed using Anaconda) command line tools, I'm trying to convert a DTED file to KML contours.
gdal_contour -f KML -i 250 e068n36.dt1 e068n36.kml

However, the KML file is wrong when loaded in Google Earth Pro: latitudes and longitudes are swapped, resulting in a wrongly located stretched contour KML near Mourmansk (it should be in middle east)

What am I doing wrong ?
KML content is consistent with the observed behavior:
  <Placemark>
    <Style><LineStyle><color>ff0000ff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle></Style>
    <ExtendedData><SchemaData schemaUrl="#contour">
        <SimpleData name="ID">0</SimpleData>
    </SchemaData></ExtendedData>
      <!-- KML expects longitude first, then latitude -->
      <!-- so everything generated below should be swapped -->
      <LineString><coordinates>37.0004166666667,68.8853787878788 37.0,68.8853787878788 36.99958333325,68.885 36.9995,68.8841666666667 37.0,68.8835416666667 37.0004166666667,68.8835416666667</coordinates></LineString>
  </Placemark>

DTED file is valid, lat/longs corner are correctly seen by gdalinfo
gdalinfo <redacted path>\e068n36.dt1
Driver: DTED/DTED Elevation Raster
Files: <redacted path>\e068n36.dt1
Size is 1201, 1201
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (67.999583333333334,37.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Point
  DTED_CompilationDate=0000
  DTED_DataEdition=00
  DTED_DigitizingSystem=
  DTED_HorizontalAccuracy=0050
  DTED_HorizontalDatum=WGS84
  DTED_MaintenanceDate=0000
  DTED_MaintenanceDescription=0000
  DTED_MatchMergeDate=0000
  DTED_MatchMergeVersion=Z
  DTED_NimaDesignator=RTEDx
  DTED_OriginLatitude=0360000N
  DTED_OriginLongitude=0680000E
  DTED_PartialCellIndicator=00
  DTED_Producer=<redacted>
  DTED_RelHorizontalAccuracy=0071
  DTED_RelVerticalAccuracy=0042
  DTED_SecurityCode_DSI=U
  DTED_SecurityCode_UHL=U
  DTED_UniqueRef_DSI=
  DTED_UniqueRef_UHL=
  DTED_VerticalAccuracy_ACC=0030
  DTED_VerticalAccuracy_UHL=0030
  DTED_VerticalDatum=E96
Corner Coordinates:
Upper Left  (  67.9995833,  37.0004167) ( 67d59'58.50"E, 37d 0' 1.50"N)
Lower Left  (  67.9995833,  35.9995833) ( 67d59'58.50"E, 35d59'58.50"N)
Upper Right (  69.0004167,  37.0004167) ( 69d 0' 1.50"E, 37d 0' 1.50"N)
Lower Right (  69.0004167,  35.9995833) ( 69d 0' 1.50"E, 35d59'58.50"N)
Center      (  68.5000000,  36.5000000) ( 68d30' 0.00"E, 36d30' 0.00"N)
Band 1 Block=1x1201 Type=Int16, ColorInterp=Undefined
  NoData Value=-32767
  Unit Type: m

Software versions:
gdal_contour --version
GDAL 3.1.3, released 2020/09/01

Google Earth Pro : 7.3.1.4507 (32-bit)


texts:
I have developed a QGIS plugin that processes an image using a bunch of other information, but it runs slowly. I would like to parallelize, as I can do with Multiprocessing.Process outside a QGIS plugin. I have seen several questions about using Multiprocessing in a QGIS plugin (for example this one), but no real answers or examples. As far as I can tell, the problem is related to management of processes, but I don't understand exactly how.

I'm using the following structure to create multiple (e.g. 5) processes:

def doprocessing():

### Test multiprocessing ###
import multiprocessing as multi
import worker

njobs = 5
jobs = []
for i in range(njobs):
    p = multi.Process(target=worker.worker)
    jobs.append(p)
    p.start()

for p in jobs:
    p.join()
    if str(p.exitcode) != '0':
        print ("Error running process_chunk as process")
        sys.exit(1)

The worker in my test just does something silly:

def worker():
"""worker function"""
#print 'Worker'

import random
import time

value = int(random.random() * 1000000)
time_of_process = time.time()

filename = "D:/Anders/SDB/multiprocessingtest_" + str(value) + ".txt"
output_file = open(filename, "w")
output_file.write(str(time_of_process))
output_file.close()

return 0

I call this from the plugin code like this (just showing you the last bit here):

    def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        import multiprocessing_test_process as process
        process.doprocessing()

When I run this through my plugin, all kinds of crazy breaks loose, with multiple new instances of QGIS starting up, and the original QGIS freezing.

If I change the first code snippet to the following and run it outside QGIS (in PyCharm), it works as it should:

def doprocessing():

### Test multiprocessing ###
import multiprocessing as multi
import worker

if __name__ == '__main__':
    njobs = 5
    jobs = []
    for i in range(njobs):
        p = multi.Process(target=worker.worker)
        jobs.append(p)
        p.start()

    for p in jobs:
        p.join()
        if str(p.exitcode) != '0':
            print ("Error running process_chunk as process")
            sys.exit(1)

doprocessing()

So it's related to the use of main or not. Any help, pointers, and preferably examples that are known to work in QGIS plugins would be great! Alternatively, if you simply know of a plugin QGIS that successfully uses Multiprocessing, please just point me to it and I will take it from there.


texts:
I have feature (polygon) was gotten by spatial query from the map.

How to change its color?

I used this code but it didn't work

IMxDocument mxDocument = application.Document as IMxDocument;
IActiveView activeView = mxDocument.ActiveView;
IRgbColor rgbColor = new RgbColorClass();
ISimpleFillSymbol fillSymbol = new SimpleFillSymbolClass();
IElement element = new PolygonElementClass();
rgbColor.Red = 255;
rgbColor.Green = 40;
rgbColor.Blue = 50;
fillSymbol.Color = rgbColor;
element.Geometry = parcel.Shape;
IFillShapeElement pFillElement = (IFillShapeElement)element;
pFillElement.Symbol = fillSymbol;
activeView.Refresh();


texts: How often do directions services update their OSM data?

texts:
Not sure if this is the right place to ask but I figured I would give it a shot. I'm working on an interactive map for a campus, I performed some edits to openstreetmaps and I can see the changes on OSM but not through third party direction services such as Mapbox and Mapzen. When I try to navigate they still follow the old paths.

I have performed smaller edits in the past and Mapbox's directions service seemed to pickup on them within a few hours to a day. I'm just wondering if there is any concrete information on this, I wasn't able to find anything when I searched.


texts: Can't change store dir or create new raster store on GeoServer 2.16. Failed to create reader from file: data/

texts:
I tryied to change raster source for existing raster store, but face it:

  Could not list layers for this store, an error occurred retrieving
  them: Failed to create reader from
  file:data/project_data/full_pop_web_2020.tif and hints Hints:
  REPOSITORY = org.geoserver.catalog.CatalogRepository@246294a2
  EXECUTOR_SERVICE =
  java.util.concurrent.ThreadPoolExecutor@70a012d[Running, pool size =
  0, active threads = 0, queued tasks = 0, completed tasks = 0] System
  defaults: STYLE_FACTORY = StyleFactoryImpl FEATURE_FACTORY =
  org.geotools.feature.LenientFeatureFactoryImpl@403c326e
  LENIENT_DATUM_SHIFT = true FILTER_FACTORY = FilterFactoryImpl
  FORCE_LONGITUDE_FIRST_AXIS_ORDER = true FORCE_AXIS_ORDER_HONORING =
  http GRID_COVERAGE_FACTORY = GridCoverageFactory TILE_ENCODING = null
  COMPARISON_TOLERANCE = 1.0E-8

I had found this post: Creating Store in GeoServer using ImageMosaic Plugin gives Error?
and I checked / tryed this:

that owner of dir/file.tif is tomcat
place a file.tif in /my_workspace/ and /data/data/ - no way
that file works well in qgis
that SRID is ok
space on server is ok

all files version were created by gdal_translate  2020_web.tif  full_web_2020.tif -co COMPRESS=LZW -co BIGTIFF=YES -co "TILED=YES" -scale 0 100 0 100 -a_srs epsg:3857

UPD full log output was added:

GeoTiffMetadata2CRSAdapter::createCoordinateSystem:Only Geographic & Projected Systems are supported.
org.geotools.data.DataSourceException: GeoTiffMetadata2CRSAdapter::createCoordinateSystem:Only Geographic & Projected Systems are supported.
    //.............................................
    ... 159 more
Caused by: java.lang.UnsupportedOperationException: GeoTiffMetadata2CRSAdapter::createCoordinateSystem:Only Geographic & Projected Systems are supported.
    at org.geotools.coverage.grid.io.imageio.geotiff.GeoTiffMetadata2CRSAdapter.createCoordinateSystem(GeoTiffMetadata2CRSAdapter.java:186)
    at org.geotools.gce.geotiff.GeoTiffReader.getHRInfo(GeoTiffReader.java:331)
    ... 160 more
2020-01-20 17:28:00,813 INFO [geoserver.web] - Getting list of coverages for saved store file:data/data_rasters/full_pop_web_2020.tif
java.lang.RuntimeException: Could not list layers for this store, an error occurred retrieving them: Failed to create reader from file:data/data_rasters/full_pop_web_2020.tif and hints Hints:
  REPOSITORY       = org.geoserver.catalog.CatalogRepository@246294a2
  EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@70a012d[Running, pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0]
System defaults:
  STYLE_FACTORY                    = StyleFactoryImpl
  FEATURE_FACTORY                  = org.geotools.feature.LenientFeatureFactoryImpl@403c326e
  LENIENT_DATUM_SHIFT              = true
  FILTER_FACTORY                   = FilterFactoryImpl
  FORCE_LONGITUDE_FIRST_AXIS_ORDER = true
  FORCE_AXIS_ORDER_HONORING        = http
  GRID_COVERAGE_FACTORY            = GridCoverageFactory
    TILE_ENCODING                  = null
  COMPARISON_TOLERANCE             = 1.0E-8

so what's the problem?


texts: Problem saving added point to GeoServer WFS-T via OpenLayers?

texts:
I want to add a point and then with WFS-T, save this to db.

I work with Openlayers and GeoServer.

saveStrategy = new OpenLayers.Strategy.Save();

wfs = new OpenLayers.Layer.Vector(
"WFS Vectore",
{
    strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
    projection: new OpenLayers.Projection("EPSG:4326"),
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        srsName: "EPSG:4326",
        url: "http://localhost:8080/geoserver/iran/wms?service=WFS",
        featureType: "population_utf-8",
        featureNS: "http://iran.kadaster.org",
        geometryName: "geom"
    })
});

I add point with two methods:

1.

wfs.addFeatures([new OpenLayers.Feature.Vector(
    new OpenLayers.Geometry.Point(
        (Math.floor(Math.random() * 360) - 180),
        (Math.floor(Math.random() * 180) - 90)
    )
)]);
saveStrategy.save();

2.

var drawFeatureController = new OpenLayers.Control.DrawFeature(wfs, OpenLayers.Handler.Point);

And then add a point with the mouse:

saveStrategy.save();

The second method works but the first method doesn't.

I need the first method.

How can I do this?

update:


texts: Export landcover area after classification - Google Earth Engine

texts:
I want to make a time series analysis of landcover change with Landsat 7 images. After I classified them using MODIS land cover classification, I want to export the change in land use. For this, I already created lists with the sum of the area per land cover type. However, I have a list for each year and do not know how to compile those together and export them.
Since by hand, this would be very time consuming.

So, how do I combine those lists for a time series analysis?

My code for this is:

// Use the MCD12 land-cover as training data.
var modisLandcover = ee.Image('MODIS/051/MCD12Q1/2011_01_01')
  .select('Land_Cover_Type_1');

// A pallete to use for visualizing landcover images.
var landcoverPalette = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  // shrub, grass, savanah
    '111149', // wetlands
    '8dc33b', // croplands
    'cc0013', // urban
    '6ca80d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
];

// A set of visualization parameters using the landcover palette.
var landcoverVisualization = {palette: landcoverPalette, min: 0, max: 17, format: 'png'};
// Center over our region of interest.
Map.centerObject(geometry, 9);
// Draw the MODIS landcover image.
Map.addLayer(modisLandcover, landcoverVisualization, 'MODIS landcover');

// Load and filter Landsat data.
var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1')
    .filterBounds(geometry)
    .filterDate('2011-01-01', '2011-12-31')
    .sort('CLOUD_COVER', false);            //sort by cloud cover ('false' indivcates descending )

// Draw the Landsat composite, visualizing true color bands.
var landsatComposite = ee.Algorithms.Landsat.simpleComposite({
  collection: l7,
  asFloat: true
});
Map.addLayer(landsatComposite, {min: 0, max: 0.3, bands: ['B3','B2','B1']}, 'Landsat composite');

// Make a training dataset by sampling the stacked images.
var training = modisLandcover.addBands(landsatComposite).sample({
  region: geometry,
  scale: 30,
  numPixels: 5000,
  seed: 0
});

// Train a classifier using the training data. USING CART() classifier
//var classifier = ee.Classifier.cart().train({
var classifier = ee.Classifier.randomForest(100).train({
  features: training,
  classProperty: 'Land_Cover_Type_1',
});

var list = ee.List
var x = 2011;

while (x < 2020) {
var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1')
    .filterBounds(geometry)
    .filterDate(x+'-01-01', x+'-12-31')
    .sort('CLOUD_COVER', false);            //sort by cloud cover ('false' indivcates descending )

// Draw the Landsat composite, visualizing true color bands.
var landsatComposite = ee.Algorithms.Landsat.simpleComposite({
  collection: l7,
  asFloat: true
});

// Apply the classifier to the original composite.
var upsampled = landsatComposite.classify(classifier);

/////////////////////////////////
// SURFC AREA CALCULATON

// get the surface covered by each classified tree
var area_classified = ee.Image.pixelArea().divide(1e6).addBands(upsampled).reduceRegion({
  reducer: ee.Reducer.sum().group(1, 'group'),
  geometry: geometry,
  scale: 30,
  maxPixels: 10e10
});
var outputReducers = ee.List(area_classified.get('groups'));
print('area of each class'+ x,outputReducers);

 x= x+1;

}


texts: Is there a way to make google charts repond to extent of map in CartoDB?

texts: I have figure out how to make a google chart responsive to a click event (view sample code here) and also responsive in an infowindow (view sample code here) but is there a way to make it respond based on the extent of the map?


texts: Update all cells in column to match value of first row?

texts:
I have a point layer; the attribute table includes x and y coordinates of original centroids associated with the points (the points were generated from the perimeters of circle features and still carry centroid data).

Since each layer contains points from two different polygons, individual point in each layer have different centroid values. I need to update the attribute table so that each row has the same value; conveniently, for each layer the desired value is from the first row (i.e. FID=1).

I'm working in ModelBuilder, so ideally this would be a model-friendly solution.

Is there an easy fix that I'm blanking on?

Possibly something in the field calculator?


texts: I have shapefiles and a tiff image that I am taking a snapshot of from a command prompt. The tiff layer always has a 3 pixel boundary. You can see the buildings and contours continue across the gap. What causes this and how can I keep if from happening?


texts:
I'm doing some tests with QField at the request of the company I work for.

I was wondering what is the best (most efficient and practical) raster format with regard to the fluidity and the loading speed while using the app on-the-field ?


texts: I have an attribute and I want to find the nearest neighbor but the issue is that the memory of the server may not be able to handle the use of all the FME transformers.  I need to take records out of a POSTGIS table and compare them to the records I have to find the nearest neighbor.  I am thinking of using ST_ClosestPoint but not sure how to write that.  If anyone has used this before it would be great to get an idea.


texts: How to update a CartoDB table at a certain time of day

texts: I have a table in CartoDB 'inland' that I would like to have automatically update each day at 8pm PST.  I know that PostgreSQL has no internal function to run an UPDATE COLUMN function at a certain time of day.  Is there a way to use Cron or PgAgent to do this?  Would I need to setup CartoDB on my local machine instead of only using the webservice?


texts: Allowing QGIS atlas to show html box in different location on some Atlas pages

texts:
I have created a QGIS Atlas, and that's contains a html box. This box is always on the same location on each atlas pages. Is there any way to custom locate the html box on some atlas pages?

I created my legend in a html box and a map view (full page) on a print composer (this means my html box is on/overlapping the mapview), the map view is controlled by atlas. In the atlas preview I have couple of pages and the html box is always on the same position on all of the pages, but sometimes it would be better to move somewhere else because it's overlapping whit some mapview features. Is it possible to do that?


texts: Need symbology that allows for symbols based on multiple fields, sometimes simultaneously (for certain features)

texts:
I've got a shapefile containing a large number of O&G parcel features.  Each feature parcel has an associated status (things like "Held by Production," "Shut-in," "Leased, Non-producing") and an associated lease expiration date.  I need to create a map that depicts this information by color-coding each parcel by its expiration date AND its status.

For example, I need to color all parcels with leases expiring in 2014 red, 2015 orange, 2016 yellow, etc.  I also need to identify those parcels with Status = "Paid Up" with diagonal lines overlaying the color fill designating expiration year.

I hope my question is clear, and thanks in advance for any feedback you can offer!

EDIT
Throughout the process of working with this problem, I've come up with another tripping point:  How does one assign symbology based only on the YEAR of a date field?  For instance, my shapefile has a DATE field with the following format:  01/01/2014.  I would like to symbolize each feature based on the YEAR only.  The symbology tool allows me to assign symbology for each individual value (1/1/2014, 1/2/2014, etc.), but not based on only the YEAR portion of the DATE field.

In fact, when I populate the symbology tool with all unique values of the DATE field and sort, the list sorts alphabetically rather than chronologically.  This leads to weirdness like 1/1/2015 appearing before 1/2/2014.


texts: How to solve "features were not reachable" warning message using v.net.distance in GRASS

texts:
I'm using GRASS 7.0.4 under Linuxmint 17.1 KDE.

I need to build a raster with the shortest distance from every pixel to a city following a road. I have two shape files, the city (point) and the road (line). I plan to compute these distances in two steps, first the straight-line distance from the pixel to the road and second the distance from there to the city following the road.

My idea is to build a point grid, compute the distances from each point of the grid to the city and transform this grid to a raster file (May be you have a better idea, could I use r.grow.distance to do this?)
So, I've built the grid with points every 500 m.

First, I've tried to solve the problem with v.distance with the to_along parameter. The problem here is that the road continues after the city so I don't know how to compute the distance up to the city and not all along the road. And I cannot clip the line representing the road because I have points to the south and to the north of the city.

After that, I've used the outputparameter of v.distance to build a vector containing the lines connecting every point in the grid to the road. I've merged this new vector with the road and obtained a network connecting all the points in the grid with the city (I think that probably here is my problem).

Then I did the following:

v.net -s input=net_clean points=grid_500 output=net1 operation=connect threshold=400 arc_layer=1 node_layer=2

v.net input=net1 points=City output=net2 operation=connect threshold=400 arc_layer=1 node_layer=3

Up to this step, everything goes fine, but when I try to compute distances:

v.net.distance input=net2 output=distances_1 from_layer=2 to_layer=3

I obtain the following message:

WARNING: 30 'from' features were not reachable

And this is doing a test with only 500 points in the grid, my original grid have almost 80,000.

I've tried with v.clean as suggested here: http://www.webrian.ch/2011/06/routing-with-openstreetmap-data-in.html and here: How to fix unreachable nodes in QGIS GRASS network?.

That reduced the original number of "not reachable" features from 117 to 30. But I still don't know how to completely solve the problem.

If you know of another tool to solve the problem, as long as it is free, I'll be glad to try it.

Any suggestion?


texts: QGIS 2.14 hangs at "Starting up the GUI" and fetching plugin repositories

texts:
All of a sudden, QGIS 2.14.4 hangs on startup on Win7 for at least 2 min. Eventually it starts without error. Hardware is not the issue (Core i7, 8GB ram, etc.). I tried reinstalling with OS4Geo installer and with standalone installer, but no change. I even cleared registry entries and .qgis2 folder. This error began after trying to connect/use a bluetooth GPS with the GPS Tools plugin. QGIS also freezes for a minute when clicking the GPS Tools button, so I suspect that this performance issue relates to the plugin in some way.

Any idea how to improve startup time?

UPDATE: Having uninstalled/reinstalled multiple times, deleted .qgis2 folder, and cleared registry, QGIS 2.14.4 still takes several minutes to start. I may try 2.16.1. Any insight into the processes running during the "Setting up the GUI" stage might be helpful.

UPDATE 2: Version 2.16 has been installed with the same result. For some reason, when my computer is connected to the network via an ethernet connection, QGIS takes 4-5min (as detailed above). However, when connected to my network via wifi, QGIS opens in about 20 seconds.


texts: Display route Events option in ArcGIS Pro

texts:
In ArcMap you can add event table and route reference feature class and when you right click on event table there is an option display route events.

I am struggling to find that option in ArcGIS Pro.
I know ArcGIS Pro has imported ArcMap functions with fancier names but it seems that this one is not included. How can I display route events? I do not want to create a actual feature class just to display route events.


texts: Listing churches with their administrative location using Overpass

texts:
I'd like to extract a list of churches from OSM with their name, their coordinates and the name of the lower administrative division they belong to. I managed to get the first ones with:

area["name"="Mayenne"]->.boundaryarea;
way(area.boundaryarea)["building"="church"];
out tags center;

I deliberately ignore nodes and relations as they don't seem relevant in this case.

I can't find a way to retrieve the name of the lower administrative division. Is there a way to do it?


texts: Generating hexagon (tessellation) fishnet inside hexagon (tesellation) fishnet using ArcGIS Pro

texts:
I'd like to create a fishnet or hexagons inside a subset area of hexagons? Basically I want a large cell size for my rural area and a small cell size for the urban area. This is the tool I've used to create the fishnet:

https://pro.arcgis.com/en/pro-app/tool-reference/data-management/generatetesellation.htm


texts: Print Composer Problems with Export as PDF

texts:
Possible Duplicate:
  QGIS Crashes when writing to PDF

My QGIS crashes during the export of maps to pdf. I already uninstalled and installed QGIS twice. The function worked for several times but quit functioning after a while. Does anyone know what this problem causes? I use a new Mac (2012). Could there be a problem with that?

Thanks a lot for your help!

Josi


texts: Aggregate Lines Based on Distance between Segments

texts:
I have a VectorTile basemap that I am optimizing, (working on reducing the size of individual tiles), the major chunk of tile size is road data (geometry and attributes). I am working on simplifying the individual road geometries based on the Zoom Level, and aggregating the roads based on zoom levels as well.
The roads are stored in a Postgis table ( individual table for each zoom level ) and I want to aggregate the roads based on the zoom level they will be displayed at. For example, In the Zoom Level 5 table I want to aggregate the roads that are within a certain distance of one another, and make a single line for the road segments:

Since it is not possible to make out individual roads at that zoom level:

How can I aggregate the road geometries based on the distance between road segments?
PS: I am working on PostGIS and QGIS on Linux, but I am open to solutions using any Open Source platform or technology


texts: Confused on gainbias.txt values using OpticalCalibration Orfeo-Toolbox

texts:
I'm currently working on a preexisting project that utilizes Orfeo Toolbox OpticalCalibration. There are different satellites that each have different data regarding solar irradiance that the project gathers from this pdf: Absolute Radiometric Calibration: 2016v0.
What I am confused about is that the pdf provides gain and offset values but it seems like the project stores the gain value through the formula absCalFactor/effectiveBandWidth and the bias value is just 0. Why didn't the project just use the values in the table?
Another extraneous question is that is the offset the same as the bias?


texts: Hiding leaflet tooltip on mobile device?

texts:
I would only display a tooltip on a device without a touchscreen. A popup and a tooltip will be displayed on your mobile by clicking on the point. I would like to see only popup

 var ubytovaci = L.geoJSON(null, {
    onEachFeature:  forEachFeature2,

    pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {icon: iubyt}).bindTooltip(feature.properties.name);
    },
    filter: function (feature, layer) {
        return (feature.properties.category == "6");
    },
});


texts: Layer adds and removes itself within a second in arcpy

texts:
I am trying to create a layer with a feature class and adding to arcmap. using below code.

import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True

inputLayer = arcpy.GetParameterAsText(0)

# get the map document
mxd = arcpy.mapping.MapDocument("CURRENT")

# get the data frame
df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

arcpy.MakeFeatureLayer_management(inputFCpath, layerpathName)

# create a new layer
newlayer = arcpy.mapping.Layer(layerpathName)

# add the layer to the map at the bottom of the TOC in data frame 0
arcpy.mapping.AddLayer(df, newlayer,"AUTO_ARRANGE")

arcpy.RefreshTOC()

it creates a layer adds in TOC and removes itself immediately from TOC.

why TOC not holding this Layer ?


texts:
I have a very extensive GPS track from which I am trying to reduce the number of plotted points. I though that I would be able to do so by grouping by a round of the timestamp /900 ( every 15 mins ), but thats not working.

Whats the correct way to reduce the resolution of the points returned.

SELECT row_to_json(featurecollection) AS route
        FROM (
               SELECT
                 'FeatureCollection'               AS type,
                 array_to_json(array_agg(feature)) AS features
               FROM (
                      SELECT
                        'Feature'                                                             AS type,
                        lg.gid                                                                AS id,
                        ST_AsGeoJSON(ST_MakeLine(lg.start_itudes3d, lg.end_itudes3d)) :: JSON AS geometry,
                        row_to_json((SELECT property
                                     FROM (SELECT
                                       round(start_time/3600),
                                             lg.start_time,
                                             lg.end_time
                                             ) AS property))                      AS properties
                      FROM ( SELECT
                               itudes_3d::GEOMETRY                                    AS start_itudes3d,
                               MIN(timecol)                              AS start_time,
                               LEAD(itudes_3d::GEOMETRY) OVER (ORDER BY timecol ASC)  AS end_itudes3d,
                               LEAD(timecol) OVER (ORDER BY timecol ASC)                AS end_time,
                               MIN(id)                                                  AS gid
                             FROM gps_data
                             WHERE device_id = :vessel_id AND
                                   timecol >= 1499643073 AND
                                   timecol <= 1502321473 AND
                                   active = TRUE AND
                                   deleted_at IS NULL
                             GROUP BY round(timecol/900), itudes_3d, timecol
                           ) AS lg
                    ) AS feature
             ) AS featurecollection


texts: Converge origins on destination leaflet geodesic

texts:
I am using Tabletop.js to map locations on a leaflet map.  With each row I add the lat lng to a feature group and then to the map. What I would also like to do is use leaflet geodesic to draw a line between all latlng pair (origins) and a single destination (orleans). Instead what I get is a line drawn between the last origin and orleans.

Here is what I have tried...

    var Geodesic = L.geodesic([], {
        weight:3,
        opacity:1,
        steps: 50
    }).addTo(map);

    var fg = L.featureGroup().addTo(map);
    var orleans = new L.LatLng(30.053206, -89.934873);
    var origin;
    Tabletop.init({
        key: code,
        callback: function (sheet, tabletop) {
            for (var i in sheet) {
                var place = sheet[i];
                marker = L.marker([place.lat, place.lon])
                    .addTo(fg)
                    .bindPopup(place.Name)
                origin = new L.LatLng(place.lat, place.lon);
                Geodesic.setLatLngs([[origin, orleans]]);
            }
        },
        simpleSheet: true
    })


texts:
I've been out of school and working in the GIS field (mostly the field of a desk) for over five years now. The one consistency I've found when working for various companies and applying for various jobs is that there is no consistency in the way our job titles correlate to the actual tasks we accomplish on a daily basis.

Does anyone have a good list of the general tasks / job description for some of the GIS job titles:

GIS Technican
GIS Analyst
GIS Administrator
GIS Developer
GIS Applications Specialist
GIS Professional
GIS Manager
GIS / IT Technican
GIS Technical Architect

I know each title varies by region and company, but a general list would be a great resource. If I have left any out, please feel free to add.


texts: Linking Data Frames with Different Coordinate Systems

texts:
I have a problem that I am hoping might have multiple solutions.  I am working with research data that was collected on an arbitrary grid, with units in m, datum point of 1000E, 1000N. The grid is rotated 34.5 degrees from geographic north. I've created a map project with our data on it, with a custom coordinate system.  I would like to use ArcGIS basemaps for background, and accurately load GPS data (UTM NAD 83) on the same map.

I have one dataset from our surveyors that gives me a set of matching data points in both our custom coordinate system and UTM NAD 83. So far I have been able to line up my maps by

Loading the site data I want on my first data frame (A) in the custom grid.
Rotating data frame A by 34.5 degrees
Adding the surveyor points to data frame A based on the site grid coordinates
adding a new data frame (B) in NAD 83 and loading basemaps and the surveyor points based on UTM coordinates
setting both frames to the same zoom (since both have units in meters)
panning until my data points line up on both maps.This works and our site data successfully lines up with geographic basemaps and GPS data.

How do I link these two layers so that I can pan and zoom in one and have the other mirror/follow it?

I tried the solution given here: Link Dataframes in ArcGIS, using data frame A as the primary/other frame, and it didn't work: it rotates data frame B to 34.5 degrees, mis-aligning my maps, and when I pan data frame A, B still won't move. Since my coordinate systems are different, and since I have a very specific alignment point, this solution does not seem applicable.

Is there a way to create a transformation for my custom grid to project it into a normal UTM NAD 83 data frame, and avoid this problem entirely? Or to somehow geo-reference one data frame to another? I want to be able to zoom and pan without re-aligning my maps every time, and to be able to use the maps in places where there are no survey points present for alignment.


texts: changing attributes or names of shapefiles

texts:
I am trying to join a subset of two similar shapefiles, one with states from the US and the second with provinces from Canada. The datasets are available here: shapefiles from Geocommons

I noticed that in the two files the attributes are slightly different. In the US map, states are labeled as "STATE_NAME" but in the Canadian map the attributes are simply "NAME". This is a problem because I cannot join the two shapefiles into a single one. Does anyone have a work around for this?

here is my code so far:

require (raster)

#load in boundaries for plotting
state <- readOGR(dsn = '/usa_state_shapefile.shp', layer = "usa_state_shapefile")
projection(state) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
# Subset US shapefile by desired states
nestates <- c("Maine", "Vermont", "Massachusetts", "New Hampshire" ,"Connecticut",
              "Rhode Island","New York","Pennsylvania", "New Jersey",
              "Maryland", "Delaware", "Virginia", "West Virginia", "North Carolina")
state.sub <- state[as.character(state@data$STATE_NAME) %in% nestates, ]
summary(state.sub)

provinces<-readOGR (dsn = '/canadian_provinces.shp', layer = "canadian_provinces")
projection(provinces) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
canprov<- c ("Quebec", "Ontario", "Newfoundland  & Labrador", "New Brunswick", "Prince Edward Island", "Nova Scotia")
provinces.sub <- provinces[as.character(provinces@data$NAME) %in% canprov,]
summary (provinces.sub)

Im guessing that if I rename the attributes to the same thing like "NAME" then I should be able to merge the two shapefiles using some sort of rbind or cbind function.


texts: How to import a module from C:\Users\myusername\.qgis2\processing\scripts

texts:
I have a user script which is located in C:\Users\'myusername'\.qgis2\processing\scripts.

I call script "plotWLayer1" containing the line

import funcPlot

User script funcPlot is also located in the same path but when I execute the script Qgis logs the following:

No module name funcPlot

I have tried adding the path using,

import os
import sys
sys.path.append(os.getcwd())

with the same result.
What am I doing wrong? My intention is to give these scripts to fellow users within my company but if I hard code sys.append() the script would have to modified per user.

Update 1:
For the moment I am using the following code under the assumption that only the username for different users in my company would change but any suggestions would still be welcomed.

import getpass
username = getpass.getuser()
import sys
exePath='C:\\Users\\'+username+'\.qgis2\\processing\\scripts'
sys.append(exePath)


texts:
I want to be able to link landmarks from a private database to their respective OSM representations (nodes, ways or relations).

One of the approaches I have seen to get permanent-ish reliable link to something found through Nominatim is:

Find it on Nominatim
Use the osm_id to download all the tags
Make an Overpass link that uses all the tags

Is there a reliable way to do that using only Nominatim? Something like a canonical query that has a good chance of always returning the same object.

Side question: Am I right in assuming that due to the transient nature of maps all I can hope for is "a good chance of always returning the same object", i.e. there is no fool-proof method of linking to the OSM representation of physical objects.


texts:
I have a shapefile with multiple polygons of different sizes and want to split this shapefile to new ones trying to achieve a determined area, but without splitting any polygon.

For example: I want shapefiles of approximately 30ha, so i would enter with my base and would results in shapefiles with polygons sums of 25ha, 29ha, 35ha, etc. The most close i could get of the 30ha without splitting the polygons.

Any idea?