texts: GRASS r.sun (in QGIS) computes higher radiation in winter

texts:
I am using a 1m DEM layer (Poland CS92, ETRS:2180). From this DEM I've computed Aspect and Slope in QGIS through grass plugins. I did a test run for 15th day of the year and 165th day of the year. The output for 165th is lower than for 15th day of the year. Wnich is physically not possible.

So I know I am doing something wrong. But I don't know what. And I've asked a couple of people, and tried different solutions, but none of those work. Here is a list of stuff I tried doing:

include lat/long rasters when running tool
trying a different DEM (Copernicus programme EUDEM)
exported raster to a different format (from ASCII to geoTIF)
set up a GRASS database, mapset and location (still running the tool from QGIS toolbox, not the GRASS one, since my grass toolbar seems to be missing r.sun)
computed slope and aspect in arcmap instead of qgis

And at this point I am out of ideas. What could I be doing wrong?? I am just plugging in the data into a ready model, there is nothing else to be done here. Does anyone have any idea what's going on? Has anyone successfully used this model, and could maybe give me some short walkthrough tips??

Probably worth mentioning that I don't know how to use GRASS and I don't have the time to learn it (tried), so it has to work through QGIS. The learning curve for GRASS is pretty steep, I am still a GIS newbie and I have a deadline, so I can't just move this project for few months later.


texts: Controlling order KML Polygons are drawn using Google Maps API?

texts:
I am creating an application that allows users to turn on and off KML layers:

The KML layers are hosted on ArcGIS Online, and I'm bringing them into the application by creating a new KML Layer variable. The KML Layers are then added to an array variable. Finally, the layers are added to the map according to whether the checkbox is toggled on.

Here is the relevant code:

layers = [];

var parcelKML = new google.maps.KmlLayer({
                url: "https://www.arcgis.com/sharing/rest/content/items/3bdd800d59254407b4e115d56fcbc9ec/data",
                map: map,
                pane: "floatPane",
                preserveViewport: true
        });
layers.push(parcelKML);

var municKML = new google.maps.KmlLayer({
            url: "https://www.arcgis.com/sharing/rest/content/items/1dcc6f19f9c24bf689235f0488d8f90a/data",
            map: map,
            pane: "floatPane",
            preserveViewport: true
        });
layers.push(municKML);

for (var i=0; i < layers.length; i++) {
        layers[i].setMap(null);
       };
    };

//Add layers to map using toggle

//Parcels
$(window).ready(function(){
        $(".ui.toggle.parcels input").click(function(){
            checked = $(".ui.toggle.parcels input").is(":checked");
            console.log(checked);
            if(checked === true) {
                layers[0].setMap(map);
            }
            else{
                layers[0].setMap(null);
            }
        });
    });

//Municipalities
$(window).ready(function(){
        $(".ui.toggle.munic input").click(function(){
            checked = $(".ui.toggle.munic input").is(":checked");
            console.log(checked);
            if(checked===true) {
                layers[1].setMap(map);
            }
            else{
                layers[1].setMap(null);
            }
        });
    });

The layers are successfully loading to the map viewer. However, I'd like the Parcels layer to always appear on top. If I toggle this layer on and then toggle on the municipalities layer, the parcels layer is now on bottom.

How would I specify the order in which the KML layers are drawn?

I'm relatively new to the Google Maps API and JavaScript.


texts: Add Label to WMS layer Using Mapbox GL JS

texts:
When We add wms layer to map Using mapbox gl js to vector map multiple layer get called and multiple time wms layer url called thats why we get many wms layer so if we have label with wms layer  when we zoom in we get multiple label of layer. so, i want to add same label externally and allow them rotate and all features of label in mapbox gl js.
map1.addSource('wms-test-source', {
'type': 'raster',
'tiles': 'url',
'tileSize': 256
});

map1.addLayer(
{
'id': 'wms-test-layer',
'type': 'raster',
'source': 'wms-test-source',
'paint': {}
}
);


texts: Extract tract shapefile from 2006 US Census

texts:
I have downloaded the 2006 census tract files.
This version does not have the shapefiles
format but the tiger format.
I am able to load the data in python using following code:

import ogr
drv_tiger = ogr.GetDriverByName('TIGER')
abc = drv_tiger.Open('AK/')

There are roughly 18 layers in abc. The layer 6 is polygon.
Since I am new to all this so I don't understand the organization that well.
Lets say I get handle of layer 6. From there how do I extract
tract shapes and finally, how do I convert them to shape files.

2006 Technical Documentation
Location of 2006 Files on FTP

Update 2:
 Based upon new answer, I am able to load the layers in postgres using ogr2ogr and shp2pgsql tool. The sample output of the polygon layer looks as follows:

    polyid | statecu | countycu | tract | block | blkgrp | state | county
--------+---------+----------+-------+-------+--------+-------+--------
      2 |       2 |       13 |   100 |  1986 |      1 |     2 |     13
      3 |       2 |       13 |   100 |  3998 |      3 |     2 |     13
      4 |       2 |       13 |   100 |  1052 |      1 |     2 |     13
      5 |       2 |       13 |   100 |  1052 |      1 |     2 |     13
      6 |       2 |       13 |   100 |  1985 |      1 |     2 |     13
      8 |       2 |       13 |   100 |  2025 |      2 |     2 |     13
      9 |       2 |       13 |   100 |  1051 |      1 |     2 |     13
     11 |       2 |       13 |   100 |  1986 |      1 |     2 |     13
     12 |       2 |       13 |   100 |  1054 |      1 |     2 |     13
     13 |       2 |       13 |   100 |  2000 |      2 |     2 |     13
     14 |       2 |       13 |   100 |  2000 |      2 |     2 |     13

Still need to figure out how to derive the shapefile given the  polyids for each tract. It seems I need to refer to complete chain, but not sure.
Possibly, one has  to merge all the polyids for block
to form the completechain for the tract.


texts: ModifiyFeature Control with an array of vectors

texts:
I want to put all the drawn vectors in an array and make editable only those which id match with a shapeId var. Do I have to add a new control in the for loop or can I setOptions of the feature to editable?

function initControls(){
 colorPoly = document.getElementById('color_picker_'+shapeId).value;

 var shapeStyle = {
    strokeOpacity: 1,
    strokeWidth: 1,
    fillColor: colorPoly,
    fillOpacity: 0.45,
 }

 var shapeStyle2 = {
    strokeOpacity: 1,
    strokeWidth: 1,
    fillColor: colorPoly,
    fillOpacity: 0.45,
    pointRadius: 5,
 }

 var styleMap = new OpenLayers.StyleMap({
    "default": shapeStyle,
    "temporary": shapeStyle,
    "select": shapeStyle,
    "vertex": shapeStyle2
 }, {extendDefault: false});

 var newShape = new OpenLayers.Layer.Vector( "Editable",{
        styleMap: styleMap
    });

 newShape.events.register("sketchcomplete", newShape, function(){
        map.addLayer(newShape);
        newShape.id = shapeId;
        createdShapes.push(newShape);
        alert(createdShapes.length)
        toggle = 1;
    });

 editingToolBar = new OpenLayers.Control.Panel({
        displayClass: 'customEditingToolbar',
        allowDepress: true
    });

 var draw = new OpenLayers.Control.DrawFeature(
    newShape, OpenLayers.Handler.Polygon,
    {
        title: "Draw Feature",
        displayClass: "olControlDrawFeaturePolygon",
        multi: true
    });
 for (var i = 0; i < createdShapes.length; i++){
        if (createdShapes[i].id == shapeId){
            var edit = new OpenLayers.Control.ModifyFeature(createdShapes[i],
            {
                title: "Modify Feature",
                displayClass: "olControlModifyFeature",
                vertexRenderIntent: "vertex"

            });
        }
    }
 editingToolBar.addControls([edit, draw]);
 map.addControl(editingToolBar);
}

What is the right way to do it? Thank you very much for your help.


texts: How to project Google Earth KML files in ArcGIS so that they display properly?

texts: I am sure my question has been answered in the past but I would really appreciate it if someone could help me with this... I have downloaded a KML point file from Google earth and I imported it into ArcMap. I cannot figure out what to do in order to project it correctly. My dataset is set to WGS 1984 UTM zone 34N and the KML file is in WGS 1984. The difference of the file's display in ArcMap and Google earth is not that dramatic but some points fall into the sea which is bad for the anaysis I am trying to perform. I tried to project the Google earth point file in WGS 1984 Web Mercator (Auxiliary Sphere), but it still didn't display properly... Can somebody please explain how it works with Google earth and  ArcGIS projections? Thank you in advance for any response :)


texts: Renaming Feature Class within versioned Feature Dataset

texts:
I am working in ArcSDE. I have several feature datasets within an Oracle Enterprise Geodatabase. When I try to rename one of the feature classes that is within a feature dataset it "falls out" of the dataset, it is now just floating in the geodatabase on its own. And now the rest of the feature dataset is inaccessible, when I try to open it I get an error saying "Layer not found" and referencing the feature class that I tried to rename.

Am I violating some rules of ArcSDE that I am not familiar with?

Thanks.


texts:
I added an L.EasyButton to my map. Now I want it - by clicking - to open a new tab with a specific address.

I tried this but it did not work:

L.easyButton, function(btn, map){
    onclick="window.open('http://www.website.com/page')"
}).addTo(map);

Any ideas?


texts: Controlling line width in QGIS Print Composer output

texts:
I am printing a 36" x 48" map (ARCH E size) and looking to have line width consistency for my contour lines.  I want to be able to define the line width as it is printed on a 36" x 48" paper at 300 DPI.

What do the absolute values of mm or pixels refer to when describing line width?  Will they stay consistent on an exported PDF despite the map's scale or size?  Can I use these absolute values of mm or pixels or must I use map units?

I have a small reference map with a colorized elevation gradient as a reference, and would like these line widths to print the same as the line widths on the big map, despite their difference in scale.  I believe I may be able to use map units, and as a function of the maps scale & size I could define a particular 'as printed' line width for any size map.

Attached is an image of the map I am working on.


texts:
I'm working with sentinel -1 image and i'm trying to compute NRPB equation as mention here:

I have applied Frost filter and then I tried to applly this normalized difference on one image only but for some reason the values I get is 0.

This is how I have tried to calculte the normalized difference:

var nrbp=imgFrost.normalizedDifference(['VH_sum','VV_sum']).rename('NRPB');
print(nrbp);
Map.addLayer(nrbp,{},'NRBP');

and this is link for the full code with geometry:

https://code.earthengine.google.com/3eb828ed8320154560d92a62e3d53755

My end goal: to calculte the differences and get real values instead of 0

Edit:
For some reason when I write it like this it works:

var VH = imgFrost.select('VH_sum');
var VV = imgFrost.select('VV_sum');
var nrbp = VH.subtract(VV).divide(VH.add(VV)).rename('NRBP');
Map.addLayer(nrbp,{},'NRBP');


texts: dataDrivenPages.pageRow is read-only or can be modified?

texts:
In ArcMap  Help  it says pageRow is read only, but in explanation it says "The index layer fields can then be read and/or modified as necessary."

(I did not succeed to use it with setValue while getValue worked)


texts: Creating animation point using openlayers and geoserver

texts:
Before getting started on creating points that are animated like this using flash
http://adventure.nationalgeographic.com/adventure/trips/americas-best-adventures/map/
I would like to know if it is possible in anyway to create animated point using openlayers and geoserver. If it is what are the different ways? I saw this example.


texts: How to get the inverse of a raster in ArcMap 10.1

texts:
I have a raster ranging from 0 to 80,000. I want to invert this, such that 80000 becomes 0.0000125.

When I use a raster calculator as (1/x) it returns 1s where the were values and blank where there was a zero. While for zero its correct the other values are wrong.

I tried this option Here 1.0 - (float(x) / 80000.), all the values are converted well only that they are scaled. I dont want this scaling, I need the exact inverses how do i go about this?


texts: Populate a field in a network with downstream end node

texts:
I am currently working on a stormwater system with the features in a geodatabase.  I would like run a tool or a script that will go through each outfall (last downstream node) and populate a field within the catch basins (upstream nodes on connected lines).

I can use the network analyst to trace/highlight the upstream nodes but, you have to place a flag for each one and then field calculate.

Any ideas?


texts: After Pasting a Feature in Pro, Why Does Move Automatically Activate?

texts:
Since moving to ArcGIS Pro 2.5 (and still existing in 2.6) when I am editing a feature and I copy and paste it, once pasting, the modify features move command is activated. This can lead to issues of accidentally moving the feature and creates a few extra steps when I am copying and pasting multiple features.
Is there anyone to deactivate this functionality, did I accidentally enable it?
My process for copying and pasting (copying a parcel to a development dataset) is as follows:

Select the parcel
Copy using the copy function in the edit toolbar of the ribbon
Paste using the paste function in the edit toolbar of the ribbon

Once I hit paste, the modify features panel comes up and move is activated on the pasted parcel.


texts: Proj4js can't handle 90/180degree transformations from 4326

texts:
So if you try to do a transform like this

Proj4js.transform( new Proj4js.Proj( 'EPSG:4326' ),
                   new Proj4js.Proj( 'EPSG:3857' ),
                   new Proj4js.Point( -180, -90 ) );

You will get: { x: -Math.PI, y: -Math.PI / 2 }

So the issue is obvious, you can't reproject the corners of 4326 to 3857.  However, Proj4js does not notify you in anyway that a math error is occurring.

Is there a general solution to these sort of reprojection problems?

So far, I'm thinking about checking for Point / ( Math.PI / 2 ) % 0 to see if division of PI is being returned then bringing the points closer to 0 so the transformation can be resolved.


texts: What's the simplest way to address "multiprocessing.pool.RemoteTraceback" in QGIS?

texts: I am trying to do semi-automatic classification (SCP) on an image but the error multiprocessing.pool.RemoteTraceback keeps terminating the process. I believe this is a Python error but I have no knowledge on Python, I'm just a beginner in QGIS and this is my first time to try SCP. What's the simplest way I can solve this error with little to no programming involved?


texts: Rebuilding Locator with ILocatorWorkspace2?

texts:
Trying to use ArcObjects to rebuild locators in an IWorkspace, but enounter 'COMException was unhandled' error when trying to use the 'rebuildlocator' method on ILocatorWorkspace2.

Here is the code I am using:

    private static IWorkspace ws;

    IWorkspaceFactory2 wsf = new SdeWorkspaceFactoryClass();
    ws = wsf.Open(propertySet, 0);

    private static void RebuildLocators()
    {
        System.Object obj = Activator.CreateInstance(Type.GetTypeFromProgID("esriLocation.LocatorManager"));

        ILocatorManager locManager = (ILocatorManager)obj;
        ILocatorWorkspace2 locWs = (ILocatorWorkspace2)locManager.GetLocatorWorkspace(ws);

        IEnumLocatorName enumlocName = locWs.get_LocatorNames(esriLocatorQuery.esriLocator, "Address");
        ILocatorName locName = enumlocName.Next();
        //IEnumLocator enumloc = locWs.get_Locators(esriLocatorQuery.esriLocator, "Address");
        //ILocator loc = (ILocator)enumloc.Next();

        while (locName != null)
        {
            locWs.RebuildLocator(locName.Name, null);
            locName = (ILocatorName)enumlocName.Next();
        }
    }

I have tried to pass both the ILocator.Name and ILocatorName.Name, but I suspect the problem might be the ITrackCAncel.

Has anyone had success with rebuilding locators with ILocateWorkspace2?


texts: How to get real lat lon values (degrees) of line geometry in QGIS?

texts:
I have a line layer in QGIS. The project CRS is WGS84 (4326) and the layer CRS too. The goal is to have the coordinates of the line nodes in a csv - one row/field per line.
Now I tried copy/paste of the attribute table to excel to get the geometry from the hidden wkt_geom column and also tried calculating a new filed with geom_to_wkt($geometry). The thing is: I do get a list of coordinate values, but they look like this:
1025757.69916104024741799 6299395.96179542783647776

WhenI expect them to look like this:
49.144401, 9.214915

This is probably a simple conversion or settings problem, but I'm stuck trying to find out what exactly I'm doing wrong here.
UPDATE:
I tried again to get the same points from my project file and from googlemaps (which is also my background tile layer for the project so should be the same)
My project:
LineString (1025848.60394998698029667 6299398.59839912224560976, 1025839.35096344409976155 6299388.28413813747465611)

google
49.144401, 9.215357 - 49.144342, 9.215275

Also this is a brandnew test project with exactly one layer containing this one line of two nodes. I set the CRS for the project, the layer and the background google maps layer to WGS 84:

sorry for the german interface, but I think you can still see the CRS I've chosen. I really don't know what I should do differently...


texts:
I need to analyze multiple GPX files with QGIS. My workflow so far is as follows:

Add an OSM tile layer to serve as a map background.
Import the GPX file as a vector layer, select the routes and waypoints layers.
Right-click the properties for the routes layer, open properties, select QML style file for lines.
Right-click the properties for the waypoints layer, open properties, select QML style file for points.

Is there a way I can set my line and point styles as the default for newly imported GPX vector layers (either globally or within the project)?

The GPX filenames differ, thus filename-related magic is not an option.


texts: query school elements using Overpass Turbo

texts:
I'm looking for the correct way to query point/polygon elements that are

name=* but neither tagged as amenity=school and building=school tags; and
elements tagged as building=school but without a parent polygon that are tagged as amenity=school (to mark the school grounds)

I expect to see a list of nodes or ways that are tagged as building=school, but are not inside a closed way tagged as amenity=school. The purpose of the query is for Q&A of school-related elements in our region.

The following is the code I use to get the initial elements I'm interested in:

[out:json]
[timeout:25];
(
  node
    [name ~ "^.school.*$",i]["building"="school"]["amenity"!~"."]({{bbox}});
  way
[name ~ "^.school.*$",i]["building"="school"]["amenity"!~"."]({{bbox}});
);
//
// how to get parent poly for each resulting element?
//
// print results:
out body;
>;
out skel qt;


texts: Converting XYZ tile request to WMS request?

texts:
I am wondering if it is possible to convert XYZ layer request
e.g
(http://dev.openlayers.org/docs/files/OpenLayers/Layer/XYZ-js.html)

to a regular WMS request? (Basically if it is possible to convert the numbers Z/X/Y to a bounding box)

I have a working WMS on the server side that can return the needed tiles, if the convertion is possible then I would convert the Z/X/Y on the server side and pass it to my WMS service to get the tile.

I need this because I have an Arcgis generated cache and it would be nice to also support the XYZ standart for the same cache.


texts: Is there a way to save locations from a Google search query? I would want to have a points shapefile based on the results. For example: "hot spring in Laguna". Of course a geojson is also welcome or other format.


texts: Linking points with lines for display in QGIS

texts:
I have 2 point shapefiles which contain lat and longs and I need to connect these with a line for display, not analysis. I am using QGIS 3.4.
To explain further, for one individual species/animal (i.e. record) I have finding location A and finding location B in a row in csv-doc - this has been added twice to QGIS to map both lat and long points independently for a single record.
I now need to, however, link for each record point A with B with lines either across shapefiles or within one shapefile. A unique ID is attached to the records. I have thousands of individual species so cannot manually digitise the connecting lines.
This is useful but I can't find the plugin - How to join many points with lines?.
Can anyone advise?


texts: Openlayers 3 - Zoom to Selected Feature (Polygon)

texts:
I have an ol.source.ServerVector with several polygons drawn as features.  I'd like for the map to zoom to the extent of the selected polygon upon selection, but I'm hitting a wall.

I was able to find this post that helped with zooming to all of the features in the source, but not the selected feature specifically

https://stackoverflow.com/questions/23682286/zoomtoextent-openlayers-3

var extent = source.getExtent();
map.getView().fitExtent(extent,map.getSize());

That works great to zoom to all of my polygons, but I can't seem to find any way to only use the selected feature's extent.


texts: ArcMap seems to be adding up values incorrectly on spatial join

texts:
I am currently trying to do a 1-1 spatial join where my merge rule should add up the lengths of polylines. However when I do the spatial join it seems to give a very large final value.

The first image is my target feature and the second image (the orange lines) are my join feature. When calculating the lengths of the orange line segments manually I get 8917. However, when I do the spatial join, I get a value of 36899. I am not sure why there is such a large disparity.


texts: ArcMap - proportion of polygon in lon/lat grid cells

texts:
I have a polygon representing a river basin on a 1 X 1 degree lat/lon grid. I need to create a boundary using the grid cells as the boundary to best represent the outline of the basin. I would like to make a condition telling ArcMap to include the 1 X 1 cells for which the polygon accounts for more than 10% of the area of that cell. I've added an image to better explain what I am aiming to do with a rough outline in blue. Any help/advice would be great.

If the above is not possible, creating a matrix of % accounted for in the polygon would still be useful and I will manually draw

EDIT: I realise the example I have drawn is 0.5 X 0.5 cells, but I need to do this for 1 X 1 grids. With the data I have the each grid starts at half degree so one cell e.g. 2.5 to 3.5N and 31.5 to 32.5E. Hopefully this makes sense


texts:
I want to combine the following 3 commands (and eventually other ones such as functions from matching.sql right after...):

routing=> select source as s from ways where gid=10213;
  s
------
 8217
(1 row)

routing=> select target as t from ways where gid=64552;
   t
-------
 44926
(1 row)

routing=> insert into test(the_geom) select ST_UNION(the_geom) from dijkstra_sp('ways', 8217, 44726);
INSERT 0 1

So I tried many things that look like

WITH depart AS (SELECT source AS s FROM ways WHERE gid=10213),
    arrivee AS (SELECT target as t FROM ways where gid=64552)
INSERT INTO test(the_geom) SELECT ST_UNION(the_geom)
FROM dijkstra_sp('ways', depart.s, arrivee.t);

but any of my experiments worked; any idea someone, please??

Edit: Thanks to dkastl, I wrote the following, certainly not optimized but still working one-line function to get a path from coordinates :

INSERT INTO test(the_geom) SELECT ST_UNION(the_geom)
    FROM dijkstra_sp('ways',
        (SELECT source from ways
where distance(the_geom, geometryFromText('POINT(5.697 45.27433)', 4326)) =
        (select min(distance(the_geom, GeometryFromText('POINT(5.697 45.27433)', 4326))) as dist from ways where the_geom && setsrid('BOX3D(0 30, 12 60)'::box3d, 4326))),
        (SELECT source
            from ways
                where distance(the_geom, GeometryFromText('POINT(5.05169 45.0482)', 4326))
                = (select min(distance(the_geom, GeometryFromText('POINT(5.05169 45.0482)', 4326))) as dist
                        from ways
                        where the_geom && setsrid('BOX3D(0 30, 12 60)'::box3d, 4326)
                    )
         )
 );


texts: Remove specific simple markers (multiple markers for some US states that are not fully connected)

texts:
This is a follow-up of How to create a gradient color expression in data defined properties

How can I remove specific centroid simple markers? For example, Florida, Kentucky, and Rhode Island all have two markers. I would like to remove one from each state.

Any ideas?


texts:
I am new to this site, and in general, to coding for what I would like to do.

I have used this code in the past, and it worked for when I was working with a map of China, but now I would like to use it for a different region of the world, and I am stuck

I have a 3 part question:

1.) I have a map of Papua New Guinea in GIF format that I would like to convert to a raster (easy)

2) Define coordinates of raster file

3) Lay over a map of Papua New Guinea and make sure it matches the raster of Papua New Guinea that I have created.

Step 1 is done, step 2 I believe is fine (please tell me if you see an error), however, there is a bug between step 2 and 3.

I am using R, and here is the code I have written so far

library(raster)
library(maps)
library(magick)
library(rgdal)

setwd( "where I am working from ")

#upload image

Species_id<-1
Species_i <- image_read(paste("Species (", Species_id, ").gif", sep=""))

Species_i_rotated<-image_rotate(Species_i,0)
image_write(Species_i_rotated, path = "Temporary.png", format = "png")

#create raster
Converted<-raster("Temporary.png", band=1)

# cut raster to fit precisely as possible the borders of Papua new guinea

ext<-c(26,445,0, 275)

Converted.cr<-crop(Converted, ext)

plot(Converted.cr)

#### This is where step 2 starts and the issue may stem from here, however, I have followed this procedure in the past, and it worked just fine
#Define the coordinates of the borders of the rasters using the sr-org projection 7852 - Australia - modified wgs84
# borders defined manually here http://spatialreference.org/ref/sr-org/7852/

xmin(Converted.cr) <- 41368.792496

xmax(Converted.cr) <- 1715419.242596

ymin(Converted.cr) <- 8696732.363608

ymax(Converted.cr) <- 9783651.873356

# This is where my code will NOT work, and I suspect its because I used this when I was working with China, and something needs to be adjusted, Nonetheless, I get an error:

crs(Converted.cr) <- "+proj=aea +lat_1=15 +lat_2=65 +lat_0=30 +lon_0=95 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_def"

newproj<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

Converted.proj <- projectRaster(Converted.cr, crs=newproj)

Error in if (maxy == miny) { : missing value where TRUE/FALSE needed

# Here is the end goal:

plot(Converted.proj)

map(region="Papua New Guinea", add=T)

here is an example of a map I am working with


texts: Retain symbology when merging layers in QGIS

texts:
I have a QGIS project with various layers for which I have assigned icons.  When I import additional data layers and then merge these layers with the existing layers (and their symbology) a new layer is created but without the symbology.

How do I merge an imported layer, without symbology, with an existing layer that does have symbology, and get a merged layer with the original symbology?


texts: How can I define a relationship between GTFS' trips and calendar/calendar_dates tables?

texts: From the Google GTFS reference guide, the 'Service_id' field of 'trips.txt' file is referenced from the calendar.txt or calendar_dates.txt file. I want to create a database for GTFS, so how I can add a foreign-key relationship here?


texts: Clip or Intersection tool - QGIS giving empty files

texts:
I am trying to clip an Ordnance Survey Master Map (OSMM) feature based on a polygon created for the sake of clipping which does not have any attributes, but the output is an empty file.

Projection of both the files is similar.

Could you please let me know how to achieve the result either using polygon or if there is any method I can use in QGIS?


texts: I need to know if is possible to insert a label in the atlas that is different for any pages. If I add a new label, for example, a description, it will reported on every page of the atlas with the same text, but I need to insert different notes.


texts: Change the geometry of an existing layer?

texts:
I work with QGIS, v. 2.0.1. I have a point layer in which I have structures that could be represented as lines or polygons, but at the moment they are represented as a group of points. I know I can create new layers with different geometry but if I do this, I don't have the attribute table of the original point layer attached to the new line/polygon layer. And if I copy the point layer, I have the attribute table but it's not the right geometry. Can I change the geometry afterwards?

I could do a table join, right? But I don't have a common attribute and the table join is not permanent is it?

I tried the "PointstoPaths" Plugin, but it's not what I'm looking for. I also found the "Points2One" Plugin which seems perfect for my problem, but somehow it doesn't work (a shapefile is created but I can't "see" anything in the map).

That's how the map looks: for lines:

and for polygons:

I have about 800 points and the attributes show name of the place, measurements, date and much more.


texts: Determine right or left side of polyline

texts:
I'm working with addresses and making a tool to check if structures are on the correct side of the street based on the house number parity matching the road parity. The roads have directionality, so I want the tool to identify if a structure is on the right or left side of the road and then do the parity check with the house number. I've been having trouble programmatically determining the right versus left side of the road. I did find Finding left and right side of line in ArcGIS for Desktop? that suggests linear referencing tools could accomplish this, but there wasn't much specific information about how to actually code the solution, and I haven't had much success making routes and extracting information on my own.
I started coding a work around that would check if all structures on one side of the street had the same parity, but I don't think this would be as accurate as if I could determine the right/left side of the road.

Does anybody have any suggestions about how to approach this?


texts: Is Fiona's collections.OrderedDict immutable?

texts:
After opening shapefile I want to use all of its attributes (geometry column, crs, etc), but want to delete the contents of the OrderedDict in the 'properties' key. I want to append new data to it, and write it back to a shapefile.
The standard syntax for deleting keys from Python dicts like so

D = {'spam': 2, 'ham': 1, 'eggs': 3}
D.popitem()
('eggs', 3)
D
{'ham': 1, 'spam': 2}

does not apply to the fiona collections object:
How can I manipulate a fiona.collection object? Or are those immutable?


texts: Visually identifying areas where older poles predominate?

texts:
I have a dataset of 300,000 poles including a number of attributes including pole age.

What I'd like to do, ideally, is be able to distinguish areas that have more old poles.

Unfortunately I can't release the data for you to play with, but I have a few tools at my disposal to try to get information to stand out!

I have a map, below, that is simply a thematic map based on pole age. This is pretty naive, but allows exploration of what's there.

How can I visualise a fairly big point dataset like this to emphasise areas that have older poles?


texts: Arcmap Maplex labeling engine bug workaround - labeling expressions won't recognize multiple spaces

texts:
I am labeling points in ArcMap using the Maplex labeling engine using a short Python script. The labels are rotated depending on an 'angle' field.

def FindLabel ( [name], [city], [num] ):
  Buffer = "\n" + " "
  if ([city] and [num] and [name]):
     return [name] + "\n" + [city] + "\n" +  [num].replace("+", "\n")
  elif ([num] and [name]):
     return [name] + "\n" + [num].replace("+", "\n")
  elif ([city] and [name]):
     return [name] + "\n" + [num].replace("+", "\n")
  elif [name]:
     return [name] + Buffer

Labels with one line appear to have an offset (about -4 y) from the end of the point's symbol, but labels with multiple lines appear centered. I wanted to add empty lines to my labeling script so that it would appear centered, but they appear to have no effect. This article notes the bug, but provides no workaround: http://support.esri.com/en/knowledgebase/techarticles/detail/30655

Does anyone have any suggestions to center labels with one line of text?

Edit:
I didn't realize this option was enabled. If anyone else has this problem, make sure you are aware of this option.


texts: PostGIS Shapefile and DBF Loader - automatically Linked dbf file

texts:
This is my system:

Windows 7

Postgresql 9.5

PGADMIN 3

PostGIS Shapefile and DBF Loader plugin

I have imported data from a DBF file correctly with the plugin.

My problem: This DBF file changes every day. Is there any option to have the file linked with my Postgresql database?


texts: Looking for a complete list of cities by state, with area code

texts: Our company is working on moving our sales assignments from area codes to ZIP codes, and we need some kind of dataset to help move existing assignments to the new system.  Does anyone here know where we could find a complete list of cities (preferably ZIP) by state, along with the associated area code?


texts: Can you create a virtual vector from zipped Ordnance Survey Open Names?

texts:
Ordnance Survey Open Names data comes as a zipped directory containing header-less CSVs for each 10km OS Grid square. The headers are in another, one line CSV. The data contains various attributes and separate X and Y coordinates.

Is it possible to create a .vrt from this data without unzipping it (zipped size < 100MB, unzipped = GBs if you want to know why!)? The bit I'm not is possible is getting it to use the header file.

I have managed to create a virtual raster of zipped OS Terrain 50 data (zip directory of gridded .asc files) using the method shown here: https://gis.stackexchange.com/a/110018/64402 with a couple of tweaks.

EDIT: Well I went all in and it didn't work so I thought I'd start simple and build up from there:

<OGRVRTDataSource>
    <OGRVRTUnionLayer name="OS Open Names">
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:27700</LayerSRS>
        <FieldStrategy>FirstLayer</FieldStrategy>
        <OGRVRTLayer name="OS_Open_Names_Header">
            <GeometryField encoding="PointFromColumns" x="GEOMETRY_X" y="GEOMETRY_Y"/>
            <SrcDataSource relativeToVRT="1">opname\DOC\OS_Open_Names_Header.csv</SrcDataSource>
        </OGRVRTLayer>
        <OGRVRTLayer name="SO82">
            <SrcDataSource relativeToVRT="1">opname\DATA\SO82.csv</SrcDataSource>
        </OGRVRTLayer>
    </OGRVRTUnionLayer>
</OGRVRTDataSource>

This is obviously only with a header file and one data file (SO82) and neither are zipped. Now loading this into QGIS, creates a point file and the attribute table has all the right headers and there are even the same number of rows as in the data file. Unfortunately, they are all empty and no points actually render (as GEOMETRY_X and GEOMETRY_Y haven't been populated). What've I done wrong?


texts: Set zoom level for display features in OpenStreetMap

texts: I´m using OpenStreetMap to show data of a specific region in Brasil. I would like to show the county names as soon I show the map, instead of in a higher level of zoom.


texts: SAR spectral signatures for oil spill detection from different sensors

texts:
I have a rather general question about supervised classification for oil spill detection. I want to apply supervised classification on the ESA SNAP dark spot detection tool results. My question is: if I create spectral signatures from an Envisat image (containing verified oil spills) will they perform well on a Sentinel - 1 image? As a pre-processing step I would apply radiometric correction (calibration) on both images. Any suggestion is welcome.
Efi


texts: PyQGIS: gdalogr:cliprasterbymasklayer new options?

texts:
I have had used gdalogr:cliprasterbymasklayer function in PyQGIS many times before but with the latest update there seems to be a lot more extra arguments and I can't figure out how to get the algorithm to function properly anymore.

processing.alghelp('gdalogr:cliprasterbymasklayer') returns the following:

  ALGORITHM: Clip raster by mask layer
  INPUT < ParameterRaster>
  MASK < ParameterVector>
  NO_DATA < ParameterString>
  ALPHA_BAND < ParameterBoolean>
  CROP_TO_CUTLINE < ParameterBoolean>
  KEEP_RESOLUTION < ParameterBoolean>
  RTYPE < ParameterSelection>
  COMPRESS < ParameterSelection>
  JPEGCOMPRESSION < ParameterNumber>
  ZLEVEL < ParameterNumber>
  PREDICTOR < ParameterNumber>
  TILED < ParameterBoolean>
  BIGTIFF < ParameterSelection>
  TFW < ParameterBoolean>
  EXTRA < ParameterString>
  OUTPUT < OutputRaster>

  RTYPE(Output raster type)     0 - Byte    1 - Int16   2 - UInt16  3 - UInt32
    4 - Int32   5 - Float32     6 - Float64

COMPRESS(GeoTIFF options. Compression type:)   0 - NONE    1 - JPEG    2 - LZW     3 - PACKBITS    4 -
 DEFLATE

BIGTIFF(Control whether the created file is a BigTIFF or a classic TIFF)    0 -     1 - YES     2 - NO  3 - IF_NEEDED   4 - IF_SAFER

Are 16 options really needed to clip a raster? Is there a simpler option via Python? I need it for batch processing.

My best attempt got an output, but there is nothing in the file.
processing.runalg('gdalogr:cliprasterbymasklayer', rstrData, shp, 'none', False, True, True, 1, 1, 1, 1, 1, True, 0, True, "", rstrOutput+str(rstrData))

Can anyone explain to me where I'm going wrong please?

https://docs.qgis.org/2.6/en/docs/user_manual/processing_algs/gdalogr/gdal_extraction/cliprasterbymasklayer.html

Docs are out of date to reflect this new algorithm.

Last version was a lot easier to control.


texts: Using Leaflet and Web Map Tile Service (WMTS)

texts:
I'm trying to use an API for rendering tiles with a Leaflet map. A link to the API looks like this and I think that the format is WMTS:

https://api.lantmateriet.se/open/topowebb-ccby/v1/wmts/token/{your token}/?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=topowebb&ST YLE=default&TILEMATRIXSET=3006&TILEMATRIX=9&TILEROW=862&TILECOL=887& FORMAT=image%2Fpng

I'm not getting the expected result when rendering the map, the tiles does not seem to be places in the correct order. Here is a jsFiddle showing the map: http://jsfiddle.net/MalinAurora/udj8ed93/2/

I guess it has something to do with the format, but I can't figure out how to solve it.


texts: Is there a paid service for accessing a PostGIS/OSM DB

texts:
I want to render a world map with mapnik and postgis, but do not want to set up my own DB, because the needed hardware and processing time for the OSM import make it nearly impossible for me. I just want to do it once or several times, but not regularly.
Is there a service I can pay for, that allows me to access a PostGIS DB for a specific time window?

Any suggestions are highly appreciated.
Thanks!


texts: Number of lines within specified distance of thousands of polygons

texts:
I have ~5000 polygons in green (districts).

I have storms that I am interested in as red lines. Each storm has a unique ID.

Looking to add fields to the districts. The fields are 25km, 50km, 100km, 150km, 200km, 250km.

Within these fields I would like to populate the number of storms for each district. So a district may have 1 at 25km but 5 at 200km (the 25km one can go in the 200km one as well).

Would prefer ArcGIS Model Builder solution but open to ArcPY or QGIS / GDAL OGR.

I could buffer the districts as these intervals and then do a unique ID intersect could but I am looking for a more elegant solution in batch.

I realize this may be a form of Generate Near Table but I cannot quite get it to run as described above.


texts: One-liner to print first feature from GPKG

texts:
To make the findings from this GDAL issue more discoverable, I'm posting here as well.

ogr2ogr has a -limit option. I couldn't find this option for ogrinfo, while ogrinfo has e.g. -where. I would like to use the same simplicity like with the -where flag, e.g. neither specify select * nor from layer.
I use it to look at actual values of a layer. -so only gives the column names and column types, and simply ogrinfo spams the terminal if you are doing it on a file with 2.5 mio features.
I'm working inside a JupyterLab and there is no way to start any GUI applications. There are only a Terminal, iPython console and Jupyter notebooks. Besides, Ctrl+C to stop the spamming might not be carried out for quite some time.
How can I do something like ogrinfo -limit 1 (ogrinfo has no -limit option) for a GPKG file?

This is not a duplicate of One-liner to print first feature from GeoJSON


texts: Apply r.neighbours or Neighbourhood analysis (LecoS plugin for QGIS) to 3000 rasters

texts:
I need to apply r.neighbours or Neighbourhood analysis (QGIS plugin LecoS (or equivalent) to 3000 rasters and define a neighbourhood of 5 x 5 and calculate the mean.

I suppose I need to create a script with a loop that goes through the raster files and applies the same analysis settings. Unfortunately I only know how to program in R. So I'm asking if anyone can show me how to create a script that does this.

I use QGIS and I would prefer to do this using the QGIS Processing Framework. But GRASS GIS or SAGA GIS are also OK for me.

Thanks in advance

EDIT: This how I did it on Ubuntu 12.04. Thanks @Curlew.

I created GRASS Location/Mapstet and batch imported the 3000 rasters using "Common formats import (r.in.gdal)" to the mapset (I recommend deactivating the option "add to layer tree"). I used one of the rasters for setting the region (g.region).

Then I opened R, set the workspace to the folder where the rasters are stored and run this:

library(raster)
library(spgrass6)
initGRASS(gisBase = "/usr/lib/grass64", home = tempdir(), gisDbase = "/path/to/gis/data/directory", location = "Location", mapset = "Mapset", override = TRUE) #Initate the previously created Location/Mapset

ls <- lapply(list.files(".","*.tif"),FUN=raster) # Import all raster to the list "ls"

for(i in seq_along(ls)){
execGRASS("r.neighbors", input=names(ls[[i]]),size=5,output="teste",flags=c("overwrite"))
execGRASS("r.out.gdal", input="teste",output=paste(paste("/path/to/output/folder",sep="",names(ls[[i]])),sep="",".tif"))
}

This solves my problem, but I'm still curious about how the equivalent Python code for QGIS Processing looks like.


texts: Spread attribute and recursively join in standalone pyqgis app

texts:
TL;DR

Is there a pyqgis equivalent to R's tidyr::spread or reshape2::cast ? If not, how do we recursively join by attribute n layers to another without erasing previous join ?

The problem

I have a layer that contains qualitative data in a given field.
I'm trying to programmatically separate the modalities of this field and then join them to a grid layer in distinct columns.

Basically, what I'm trying to do resembles functions such as tidyr::spread or reshape2::cast for those familiar to data wrangling with R. Only I'm trying to do this in a standalone pyqgis application.

Sample

The is meant to work with any given ogr vector that has areal qualitative data. Therefore, the program must loop to get each modalities in a given data field

Data table (data.shp) after intersection with grid

id | data | area
 1 | cat1 | 123
 1 | cat. | 456
 1 | catn | Null
.. | cat1 | 012
.. | cat. | 345
.. | catn | 678
 n | cat1 | 890
 n | cat. | 123
 n | catn | 456

Expected result

id | cat1  | cat. | catn
 1 | 123   | 456  | 789
.. | 012   | 345  | 678
 n | Null  | 123  | 456

My almost but not quite working solution

The best I can come up with is this :

Get feature and append to a list with data.getFeatures()
Split data layer with processing.runalg("qgis:splitvectorlayer")
This creates n shapefiles in a temporary folder that I destroy at the end of the program.
In a for loop on point 1. list, load/Edit the splitted files by reconstructing the path with .format, storing each QgsVectorLayer() in another list and finally editing the field name with lyr.startEditing() lyr.renameAttribute() and lyr.commitChanges()

This all works fine.

Final step is to recursively join splitted layers onto the grid layer. The following doesn't work :
In a for loop on point 3. list (modalLyr), I reconstruct the field name and attempt an addJoin(), like this :

```

for lyr in modalLyr:
areaField = "sup_{}".format(modalLyr.index(lyr))

joinObject = QgsVectorJoinInfo()
joinObject.joinLayerId = lyr.id()
joinObject.joinFieldName = "id"
joinObject.targetFieldName = "id"
joinObject.setJoinFieldNamesSubset([areaField])
joinObject.memoryCache = True
grid.addJoin(joinObject)

```

The problem with this is that I only get one additional field in my output, instead of n areaFields...


texts:
I am trying to add a point with GeoAlchemy 2. However I am getting getting error on insert action saying

---sqlalchemy.exc.InternalError: (psycopg2.InternalError) parse error - invalid geometry
HINT:  "POINT(39.9130572,3" <-- parse error at position 18 within geometry
 [SQL: 'INSERT INTO spots (name, geom) VALUES (%(name)s, ST_GeomFromEWKT(%(geom)s)) RETURNING spots.id'] [parameters: {'name': 'testspot', 'geom': 'POINT(39.9130572,32.7892743)'}]---

I check the spots table and its geometry appears as point.Could anyone tell me where is the problem ?

class Spot(Base):
    tablename = 'spots'
    id = Column(Integer, primary_key=True)
    name = Column(String)
    geom = Column(Geometry('POINT'))

Spot.table.create(engine)

if not session.query(Spot).filter_by(name='testspot').first():
    spot = Spot(name='testspot', geom='POINT(39.9130572,32.7892743)')
    print(spot.geom)
    session.add(spot)
    session.commit()
else:
    spot=session.query(Spot).filter_by(name='testspot').first()


texts: How to combine raster and attribute table in QGIS?

texts:
I have some rasters representing string categorical variables. Each raster contains integers, used to translate the categorical  variables. I have also the .csv tables giving the meanings of the integers used in the rasters (1='Corn'; 2='Soybean', etc for example).

Is there a way to combine the rasters and the tables in order to:

a) Display the strings and not the integers when I use the "Identify Features" tool?

b) Use the strings instead of the integers in the Raster Calculator tool?

I am a QGIS 2.4 user on Windows.


texts: Algorithm to find low/high atmospheric pressure systems in gridded raster data

texts:
What's the commonly used algorithm for finding center positions and values of low/high atmospheric pressure systems in gridded raster data (GRIB/NetCDF/GeoTIFF)?
First I tried a naive approach to iterate through the data array, and find points with smaller/larger value than all their neighbors. This gives surprisingly good results, however it also gives lots of false positives, especially near the equator. I need to filter the results somehow, but I'm out of ideas how. It seems there needs to be some tolerance setting.

I tried blurring the data, it helps, but still it is not perfect. Also, it often shifts the positions outside of the smallest contour, so it can be perceived as incorrect.
I already use marching squares to find contours. A possible approach could be to chain the functions together, find min/max in every smallest closed contour.
However, I still feel that it should be possible somehow to implement separately from contours. Is there a way? Ideally in JS, but I welcome suggestions in any language or tool if it mentions which algorithm does it use, so that I can reimplement it in JS.


texts:
We have ArcMap 10.6 at work, and would like create choropleth maps to visualize the  the change in profit (%) that adjusting different business variables would have.

**Note that the actual change in profit probably range from -30% to 40%, but I'd still like to have bands like the ones below that range from -100% to 100%, so that we can easily compare different maps with different scenarios.

For the maps, we'd like to have set 10% profit bands ranging from -100% to 100%.

Something like:

x<-100%
-100%<=x<-90%
-90%<=x<-80%
...
80%<=x<90%
90%<=x<100%
x>100%

How can I accomplish this? I have right-clicked the layer, gone to properties, symbology, picked graduated colors, and picked the variable. When I click classify to create my bands, I can't seem to create the bands that I want.

Also, if I were able to create the bands that I want, is there a way to quickly import them to another map, to avoid rework in getting the bands the way I want them for a new map?


texts: vary bubble size and color in CartoDB Torque Cat

texts:
I'm trying to make an animated map using the Torque Cat feature in CartoDB editor.

But I could not figure out how to edit the cartoCSS to vary both bubble size and color based on my dataset. I want to vary the bubble size based on a column with numbers (number of A), and also vary the color based on another column with strings (country name).

According to the documentation, it is not possible to use more than one column.

  cluster does not contain values for all the columns, you can only use
  value variable you would need to use an aggregation function for
  strings

A blog post by CartoDB's Andrew Hill on Torque examples did not address such issue too.

Anyone has experience working on this issue?


texts: Specifying different margins on each side of map in a QGIS report

texts:
When generating a report with QGIS it is possible to define a margin around the selected atlas feature when drawing a map. Unfortunately it isn't possible to define different margins for top, bottom, left and right. I am certain that there must be a way to do this with a python function, but I am utterly lost and my googling has not come up with anyone in a similar boat.
The reason I want to do this, is that I have a bunch of other stuff with semi-transparent backgrounds arranged on top of the map, and I want these to be shown on top of a background raster layer from the map, but not on top of the selected feature. So ideally I would define the margins in page units rather than map units.


texts: Pairwise shortest distance between points in a road network

texts:
1. Problem
I have the coordinates for over 5000 points (cities) that lie over a road shape file, as partially illustrated below:

I'm working on R, and I need to find the shortest distance between all pairs of points along the roads.
2. What I've already tried, without success:

stplanr package: had a problem converting the full road shape file to a SpatialLinesNetwork object ("non finite line lengths");
osrm package: can't compute all distance pairs (over 30 million) due to server restrictions (up to 10k pairs).

3. Sample data
I provide a sample data here: https://github.com/cassianord/example_data/blob/main/example_data.7z?raw=true


texts: Batch process to get single band image from stacked image using ArcGIS, ERDAS or ENVI?

texts: Can you help me to extract a single band image from stacked multiband image. Totally, 260 image has been stacked. When I do to this, it will take a long time. So I'd like to know there is any batch process to extract single band image in ArcGIS, ERDAS or ENVI.


texts: Changing database connection parameters in QGIS projects stored in PostgreSQL

texts:
We have stored our QGIS projects in a PostgreSQL database.

Now we have to move the server that hosts the database that many projects access. Is there a way to change the connection parameters within each project afterwards?

For projects stored outside of the database, the XML can be easily changed in a text editor, but unfortunately the projects stored in the database are in binary format.


texts: Saving to different schema using Osmosis?

texts:
Does anyone know how to use osmosis to save to a different schema?

I used the pgsnapshot_schema_0.6.sql to create a blank data table, but it always creates it in the public schema.  I need to create these table in the 'gis' schema.

Can this be done?


texts: .NET Components triggering ITools on ArcEngine/Objects

texts:
I'm building a simple map interface on a .net component that I need to redistribute around.

Well, my initial thoughts was to create the tools just the way I'm used to do it in ArcMap, inheriting BaseTool and then set the MapControl.CurrentTool = tool in my onclick form events.

What is the best way to accomplish this? This one?

Also, how can I pass the hook to those tools? I'm using this class:

public class CommandPoolManager
{
    public ICommandPool _CommandPool;

    public ICommandPoolEdit GetCommandPoolEdit
    {
        get { return (ICommandPoolEdit)_CommandPool; }
    }

    public ICommandPool GetCommandPool
    {
        get { return this._CommandPool; }
    }

    public CommandPoolManager(object hook)
    {
        this._CommandPool = new CommandPoolClass();
        ICommandPoolEdit commandPoolEdit = (ICommandPoolEdit)this._CommandPool;
        commandPoolEdit.SetHook(hook);
    }

    public void AddCommandToPool(ICommand command)
    {
        ICommandPoolEdit commandPoolEdit = GetCommandPoolEdit;
        commandPoolEdit.AddCommand(command,null);
        commandPoolEdit.CallOnCreate(command);
    }

    public void AddCommandToPool(ICommand command,UID uid)
    {
        ICommandPoolEdit commandPoolEdit = GetCommandPoolEdit;
        commandPoolEdit.AddCommand(command, uid);
        commandPoolEdit.CallOnCreate(command);
    }
}

Would that be the correct way to use CommandPool? I'm without options for testing this, since it's a component to be distributed. The real question is: Am I on the right track?


texts:
I'm looking for point pattern analysis algorithms and literature on point pattern analysis, online resources and book titles welcome. The topics may range from general descriptions of algorithms to concrete use cases in any field of research.

Update Jul 31, 15:54:

I'd be especially interested in local methods to detect linear features in point patterns. e.g. create road/track geometries from gps points.


texts:
How can I convert a point in a given CRS I know, to ECEF, using Proj4?

Is there an EPSG code for ECEF coordinates? Or is there a way to create a proj4 transform object to convert from and to ECEF coordinates?


texts: Pure Python Library for Geometry Operations

texts: Does anyone know of a pure Python (no C or Java dependencies) library for doing geometry operations such as buffer, union, intersects, etc?  I'm looking for something that is Google AppEngine compatible.


texts: Mosaic with LISS-III RESOURCESAT-1 images

texts:
I'm having problems with dark edges in normal or virtual catalog (.vrt) with LISIII LISS-III RESOURCESAT-1 mosaic. If I put no data value 0 in properties, change border by white (see figure) and I want a superposition between satellite images in my final geoTif image. Can someone help me?

My code:

gdalbuildvrt "C:/Users/ASANTOS/Dropbox/GIS_formiga_Arnhold/Data analysis/test.vrt" "C:/Users/ASANTOS/Dropbox/GIS_formiga_Arnhold/Data analysis/P6_LIS3_20120311_324_100_L2_BAND2_FINAL.tif" "C:/Users/ASANTOS/Dropbox/GIS_formiga_Arnhold/Data analysis/P6_LIS3_20120316_325_100_L2_BAND2_FINAL.tif"


texts: Multipanel with raster (spplot) and scatter plot in R

texts:
I have one raster and one scatter plot that I want to put together in a single multipanel plot. However, the par(mfrow=c()) that usually works to put together different normal plots, it doesn't work when I am trying to visualize a raster with the sp::spplot function. I am looking for a way to put both types of plots together in a single figure

This an example code what I have been trying to do, but end it ends with the plot in separate figures.

library(sp)

r <- raster(system.file("external/test.grd", package="raster"))
x <- rnorm(100, 1, 2)
y <- rnorm(100, 1, 5)

par(mfrow=c(1,2))
spplot(r)
plot(x, y)


texts: Random behaviour with rule base labels (v3.4.1)

texts:
I have a layer with rule based labelling enabled. There are two label styles: a base style for all labels (no filter) and a highlight style which is intended to apply to features with a particular attribute value only ("Class" = 'Review').

Sometimes it works, sometimes it doesn't. Each time the map canvas is refreshed each feature with the 'Review' value randomly chooses either the base label style or the intended highlight label style.

Have I failed to set this up properly or is this a bug?

Rule based label implementation with style applied randomly to one feature.

Refresh and this time the rule failed to apply to anything.

Another refresh and this time the rule randomly applied to another feature.


texts: Convert WGS 84 Coordinates to Arc 1950 coordinates


texts:
I'm trying to make an alert window show the current mouse coordinates in openlayers 2 using the google maps projection EPSG:4326. Essentially like this post but in openlayers2.

But its not working out like that. I get . . .

Heres the code inside the click event.

var point = map.getLonLatFromPixel(e.xy);
                    myXY = new OpenLayers.LonLat(point.lon,point.lat).transform('EPSG:4326', 'EPSG:102113');
                    alert(myXY);


texts: arcpy.SelectLayerByAttribute not working with an integer saved in a variable

texts:
I have a long code written and this is the last little bit I need and for whichever reason it does not work. This is the simple code.

k=2
arcpy.SelectLayerByAttribute_management("poligon", "NEW_SELECTION", "Id=k")

I tried a lot of different apostrphe, brackets and everything, nothing works. It says a column was specified that does not exist, Error 000358. But when I replace k with a number everything is OK. I'm at a loss of what is wrong.


texts: Displaying numerical attributes as strings in ArcGIS Online?

texts:
I have several columns of research data that are coded 0-2. For example:

0=No monkeys in state
1=Some monkeys in state
2=Only monkeys in state

When I actually bring this map up to ArcGIS Online, obviously I don't want someone to click on a state and get a popup telling them just an integer. I would like some way of writing a conditional statement that will display whatever string I want instead its numerical code.


texts: Getting latitudes and longitudes in between start and end coordinates and spaced at equal distance

texts:
I have a csv file which has the latitude and longitudes values , I need to find the coordinates in between the start and the end latitude and longitude which must be equally spaced at the interval of 10 meters.(to refer here)

The csv file is as follows:

 id latitude longitude
  0   2.456     4.566
  1   3.457     6.455
  2   4.890     7.787
  3   5.679     9.852

From the given csv the start lat and long values are of id 0 and end coordinates are of id 1 .then for second time the start coordinates must be of id 2 and the end point must be of id 3 and so on. These must iteratively take the values from the csv and find the in-between coordinates from the start and end latitude and longitude values.

The output must be in this format:

 latitude longitude in_between_lat in_between_lon
   2.456     4.566     2.456          4.566
                       2.466          4.576
                       2.476          4.586
                       2.486          4.596
   3.457     6.455     3.457          6.455
   4.809     7.787     4.809          7.787
                       4.819          7.797
                       4.829          7.807
   5.679     9.852     5.679          9.852

From these we can conclude  the in-between coordinates between start point 2.456, 4.566 and end points 3.457,6.455 and so on (here the values taken for the in-between lat and long are random).


texts: Failed connection to Earth Engine servers

texts:
I'm trying to export a multi-layer image to drive on Earth Engine, but the task fails with the following error:

  Error creating or submitting task. Failed to contact Earth Engine
  servers. Please check your connection, firewall, or brows

and then the message cuts off.

I've tried getting around this by logging out and logging back in, restarting my computer, and even using separate Earth Engine accounts both on my home wifi and university VPN for internet, so I don't think it's specific to one Earth Engine account or internet situation. Also tried over two days now, so I don't think it's a temporary server problem either.

Repeatable code below; wondering if anyone can provide insight.

// Generate Region of Interest
var ROI = ee.FeatureCollection('USDOS/LSIB/2013')
  .filter(ee.Filter.eq('cc','GL'));

// Compile the data
var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2005-01-01', '2005-12-31'));

// Select the data you're interested in
var ndvi = dataset.select('NDVI');

// Prep data for export
// https://gis.stackexchange.com/a/254778/67264
var stackCollection = function(collection) {
  var first = ee.Image(collection.first()).select([]);
  var appendBands = function(image, previous) {
      var dateString = ee.Date(image.get('system:time_start')).format('yyyy-MM-dd');
      return ee.Image(previous).addBands(image.rename(dateString));
  };
  return ee.Image(collection.iterate(appendBands, first));
};

var ndvi_img = stackCollection(ndvi);
print("NDVI Image", ndvi_img);

// Export the data
// Export a cloud-optimized GeoTIFF.
// See https://developers.google.com/earth-engine/exporting
Export.image.toDrive({
  image: ndvi_img,
  description: 'MOD13Q1_NDVI',
  scale: 250,
  region: ROI,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});


texts:
Currently I am adding a basemap from stamen.com to my plots as it is explained here: http://geopandas.org/gallery/plotting_basemap_background.html .

Is it also possible to load a WMS from another source as basemap?

I'm interested in this one: http://sgx.geodatenzentrum.de/wms_topplus_open


texts: Reducing the size/extents of a QGIS project created from geodatabase (.gdb)?

texts:
I have imported a .gdb geodatabase into QGIS by dragging and dropping the .gdb file into QGIS layers region. The project covers a large amount of territory and there are a lot of layers.
Due to the project's extents and size, QGIS is very glitchy and slow to the point of almost being unusable.
Is there a way to zoom into a smaller region and trim/cut the project extents to this area? Perhaps exporting it to a new QGIS project file that can be worked on vs. the large project?


texts: is there any examples or use cases of using GIS in outer space and in no-earth related use ?


texts: Making data table using linear kriging interpolation with variogram in R

texts:
I have a table of data with such structure:

longitude_for_data = c(32,68,89,145,176, -14, -42)
latitude_for_data  = c(22,98,21,13 , 16,-134,-102)
Z_for_data         = c(30,40,60,70 , 20,  40, 100)
data_to_analise    = cbind(longitude_for_data,latitude_for_data,Z_for_data)

I want to make grid and fill it using linear kriging with variogram(as a variant -simple kriging):

longitude    = seq(-180,180,by=5.0)
latitude     = seq(-180,180,by=2.5)
lat_long_grid= expand.grid(longitude,latitude)

as a result I want to get something like that:

   lon  lat Z-kriging
1 -180 -180 32
2 -175 -180 35
3 -170 -180 38
4 -165 -180 39
5 -160 -180 40
6 -155 -180 41

How should I use kriging function to recieve that data table?


texts: Why does QGIS only intersect one part of the layer?

texts:
I have a QGIS layer with information about administrative boundaries at sub-municipality level. I want to extract only land area in each sub-municipality using information from another layer. I've followed the suggestions here: How to create a new layer from overlap between two layers?.

Each row in the "administrative division" layer attribute table corresponds to one sub-municipality, and each row in the "land area" layer attribute table corresponds to one municipality (for some reason municipality data covers only land area, while sub-municipality data also includes water).

When I use the intersect tool, the created new layer only covers some parts of the layer (i.e., some municipalities, as displayed in picture below*). The pink/purple part is exactly the result I want, but I want it for the entire layer. What is the problem? I've ensured that all rows are selected in both layers.

This is my first time working with GIS-data, so forgive me if it's a stupid question.

**Data privacy note: I've edited the picture so these are not the true administrative boundaries.*


texts: How to register dll extension with EsriRegAsm.exe?

texts:
I am having some issues with registering a dll:

run DOS as Administrator
CD "C:\Program Files (x86)\Common Files\ArcGIS\bin\" (where the ESRIRegAsm.exe is)
I attempt to register the dll as follows:

Syntax:  ESRIRegAsm.exe assembly_name /p:product_name

ESRIRegAsm.exe "C:\Program Files (x86)\MyExtension\myextension.dll" /p:desktop

I get a warning raised by the registration utility that says: Specified path not found.

However, the path is positively correct.

What am I doing wrong?


texts: CadTools plugin functionality in QGIS 3 beyond Advanced Digitising Panel?

texts:
I've recently migrated to QGIS 3.6.0, and I've realised that the wonderful CadTools plugin is not compatible with QGIS 3.

I tried to use the Advanced Digitising Panel, but it's cumbersome and not practical at all, especially when it comes to drawing land-use plans. For instance, with CadTools, I could easily create arc intersections, but with the Advanced Digitising Panel, I couldn't find a way to do it without losing my mind. The CadTools plugin is super efficient and straightforward to use, and I already miss it so much. I resorted to having both QGIS 2.18.28 and 3.6.0 installed on my computer as a workaround...

Does anyone know if there's a replacement for this plugin or if it'll be converted for QGIS 3?

I contacted the plugin author, but I haven't received any answer, and the plugin's Github seems to be dead.


texts: Is there a standard way to apply different field values to linestring vertices?

texts:
In QGIS I have several 3D points that I am creating a linestring from and each point could have a different value in one of the fields like "type". I want to somehow create the linestring and be able to get individual info for each GPS point.    The 'type' only has 4 choices A,B,C or D.  The line could have 50 vertices.   I looked at multipoint also and it also seems to allow only the same values in each field for each point in the multipoint feature.

I'll have hundreds of lines when this is finished.

Is there another data type or a way to group points into a collection from which I can make a line?


texts: Calculate azimuth of vector for QGIS 2.0 Python plugin

texts:
I calculate the azimuth of a vector and update a field in my QGIS 2.0 Python 2.7 plugin. This works currently, but it takes much longer than i'd like to. I believe it is because i'm querying the data source on every call. Any suggestions on updating the fields all in one go as was suggested in another one of my questions? I'm having issues trying to adapt the suggestion for this particular evaluate function.

    self.rotateLayer.startEditing()
    bear_idx = self.rotateLayer.fieldNameIndex('angleRot')
    bearCalc = QgsExpression('(atan((xat(-1)-xat(0))/(yat(-1)-yat(0)))) * 180/3.14159 + (180 *(((yat(-1)-yat(0)) < 0) + (((xat(-1)-xat(0)) < 0 AND (yat(-1) - yat(0)) >0)*2)))')
    bearCalc.prepare(self.rotateLayer.pendingFields())
    for fields in self.rotateLayer.getFeatures():
        fields[bear_idx] = bearCalc.evaluate(fields)
        self.rotateLayer.updateFeature(fields)
    self.rotateLayer.commitChanges()


texts: I have calculated a temperature map from Landsat Image. This map is correct but I the cells can be seen obviously. I have changed the cell size from 30m to match my other products which are 1m. I used resampling tool and the cell size is indeed changed when I check the source data but not in visualization. How can I make these apparent cell borders disappear?


texts: DTED/GeoTiff layer: Get elevation under mouse pointer

texts:
I have a simple desktop app with a map control. The map data that will be available to show will be either DTED files (Military Analyst's catalog / MADTEDLayer), or GeoTiffs (created by gdal). I want the user to be able to continuously see the elevation under the position of the mouse cursor. Is it possible to extract this information from the map layer somehow?

PS. I suppose that GeoTiffs have some kind of elevation data embedded since I opened them with GlobalMapper and saw height information under cursor.


texts: Using field names on output filename in arcpy

texts:
I am trying to write a loop in arcpy that will first make a selection by attribute in one feature class and then make a selection by location on another. Then I want to copy the selected features into a gdb and move to the next row in the table.

This works just fine when I use a numeric uniqueID, but I want to have the output filename retrieve a value from the attribute table (A string with the feature's name) instead of giving me "Pand + string_x" which is just a number.

This is the script:

x=1

while x < 389:

string_x = str(x)

x_input = "UniqueID = " + string_x

x_output = r"C:\temp\BaseMap\working_data\BAG.gdb\Pand_" + string_x

arcpy.SelectLayerByAttribute_management("gemeente2017","NEW_SELECTION", x_input)

arcpy.SelectLayerByLocation_management("Pand","HAVE_THEIR_CENTER_IN","gemeente2017","","NEW_SELECTION","NOT_INVERT")

arcpy.CopyFeatures_management("Pand",x_output, "")

x += 1

I'm guessing the solution will have something to do with Field delimiters, but I'm quite new at this and have no idea how to incorporate that into my code.
Using ArcGIS 10.5 and Python 2.7.12.


texts: OTF reprojection and reprojecting layer not working in QGIS

texts:
more troubles in QGIS.
I georeferenced this raster in Corrego Alegre UTM23S and loaded up, so project CRS automatically changed to that. So far so good. I need it to be however in WGS84 UTM 23S. None of the two possible ways will do. Enabling OTF reprojection and setting project CRS to WGS84UTM makes the raster layer disappear. Even after "zoom to layer". It's just not visible.
The other way would be to save the layer in a different CRS with right-click-save as. Doesn't work either. I just did as other times before, but the file is not saved. I open the folder where I sent it to, and it's not there.
Please help.


texts:
As far as I understand WPS (web processing service) can be used to execute geospatial operations on the server side.

Is it possible to create scripts with WPS?

For example something like:
'check in a directory if there is a new file' - if so , 'process this file', 'store the result in the database'

I assume that I'm expecting too much, but what is the best solution for an script when using qgis server or geoserver?

D3


texts: Viewshed analysis without observer point?

texts:
I want to find the most visible point within an area..
I don't have an observer point, I want to find it.

based on the highest points, I've created scenarios and compare..
However there might be an automatic way of doing it?

I have access to ArcGIS desktop (spatial analyst), QGIS or Global mapper..
the software package doesn't really matter.. I just want to know if is there something already there.


texts: Filling in patterns into attribute table column of ArcGIS?

texts:
Attribute tables in ArcGIS are not nearly as easily editable in comparison to something like Excel.

For example, you have sorted the table based on some column value, and now want to fill in the column beside it with 2,4,6,8... etc. descending.  In Excel this is extremely easy - just start the pattern, highlight it and drag down. Done.

If I wanted to do something like this natively in ArcGIS, from what I know of ArcGIS, python scripting is required.

Of course I could also look at the table through Access/Excel/some spreadsheet, and do the operation from there much easier, save it, and open it back up in ArcGIS.

However, ESRI does not officially recommend such tricks.

Really, it would just be great if attribute tables in ArcGIS (and all GIS for that matter) had the same type of editability as Excel, but I'm sure there are many spatial database issues with even beginning to think about implementing such a concept.

So is the "best practice" for filling in patterns into tables really using python scripting?

I can write python scripts, but have a harder time managing what I've done in the past (maybe I just need to get better at storing "common" scripts?).

Or is there officially recommended ways that avoid scripting by using tools/extensions in ArcGIS, or some other software?


texts:
I was wondering initially how one can rename a user profile in QGIS (3.4). It seems there is no feature for that in the GUI.

Eventually I wondered if it's possible to duplicate a user profile. I did not see any topic about this anywhere on the web but I concluded it might simply take to copy-paste the folder that contains the user profile data.


texts: Calculate statistics for certain areas of a raster

texts: I want to calculate the value of certain areas of my raster, which are determined by a polygon layer. I can't use the zonal statistics tool, I receive an error about inputs. How should I do it?(the green spots is my polygon layer)


texts: Lines To Polygons and Polygonize attribute table problem

texts:
On my systems, QGIS 2.18 64-bit, Windows 10, attribute tables produced by the Lines To Polygons tool and Polygonize tool differ from the table of the source layer.

The source table has two fields:

ID--type Integer64, length 18;

PARCEL--type String, length 32.

All fields in the output table are type string, length 255.

Is there some setting or step I am missing to produce output tables that have the same field type and length as the source table?

The Lines to Polygons tool states: "The attribute table of the output layer is the same as the one of the input line layer."

I have tested on several other shapefiles.

Also, I get similar results with the Refactor Fields tool, i.e., all output table features are type string, length 255, which I addressed in this thread QGIS Refactor Fields unwanted behavior.

Update 10/8/17
Appears to be bug in OGIS 2.18. Upgraded to latest version 2.18.13 today and tested 'Lines to Polygons' tool with same problem results.

Then downloaded and installed QGIS 2.14.19 and ran same test. Results are as desired and expected: output layer attribute table is same as input layer.

Update 10/10/17
Today I filed Bug report #17266 https://issues.qgis.org/issues/17266.


texts: Combining sections of a polyline in QGIS whilst summing a attribute table column

texts:
ive got to condense a sewer layer from junctions (small sections of pipe between manholes)to one section of pipe between the manholes and sum the lengths of all the sub sections.

Does anyone know of a way to do this in QGIS?


texts: Specify field of raster in Raster Calculator expression of ArcMap

texts:
In Raster Calculator I need to apply the next expression:
Con("raster1%","field3")==Integer, Lookup("%raster2%","field3"),Lookup("raster1%","field3"))

or more specifically:
Con("%alteracionMapeada_raster%"==8,Lookup("%swir%","CLASS_NAME"),Lookup("%alteracionMapeada_raster%","OCURREN16"))

The "%raster2%","field3" and "raster1%","field3" are categorical values so I need to get those values (form "%raster2%","field3" and "raster1%","field3") in the output raster cells. The problem is that the results is just one numerical field (no other fields with categorical).
I have already checked Specify a field in attribute table with raster calculator and the answer about the lookup tool has very promising but didn't work.


texts: Font size in text annotation when exporting QGIS layouts to PDF or JPEG?

texts:
When I try to export a map from the Layout Manager in QGIS, all the elements apart from the text annotation come through fine. These come through as tiny objects, with tiny font.

I could put labels on in the layout manager, rather than on the project window, but as far as I can see, there is no option for using the handy call-out boxes that are available in the main project window.


texts:
I'd like to make a 3D population density map about a country.

I've got all of the population data of the towns and villages etc. and I've got shp about these areas.

I'd like to make a map like this:

I tried to do this with QGIS ThreeJS plugin. This is the result:

I'd like to add state vector layer above the 3D layer. Is this possible with that plugin?

Can somebody tell me what tool could I use to do this?


texts: Searching for GIS data sources relating to the Soviet GULAG system

texts: I'm desperately searching for GIS data related to the Soviet prison camp system, the GULAG. There are quite a few very good maps out there, but I cannot seem to find data in order to make my own. Thinking of locations of the camps, administrative divisions of the penal system etc. Any ideas, anyone?


texts: Which DEM format is efficient .tiff or .hgt

texts:
I know this is a vague question, but I've no expertise in GIS world. I was building this terrain map for a game like simulation project that I'm invested in. Data accuracy is least or no concern, but the quality of the terrain it generated. I'm not storing 27GB of data, instead I'm procedurely generating the height maps, from low quality SRTM30 (couple of GB) as I need them. I use mainly GDAL and interpolation (Kriging), numpy and some Scikit-learn. But I was unable to do this in real time.

Which DEM format is most suitable for GDAL or an expensive loop to work with, to produce close real time result.

Are there other choices?


texts: Looking for location intelligence tools and user experience reports

texts:
I am looking for location intelligence tools (= combination of business intelligence + GIS). I had a look at several options: SpagoBI, Microstrategy + VisualCrossing, Yellowfin, and Map Intelligence.
Can anyone tell me about their experience using these tools? Or maybe I missed some other options?


texts:
I am using ArcGIS Runtime SDK for Android and it supports wmslayer but there is no sample or code for selection feature with using wmslayer. My app shows the map getting data from my wms server.

Do you know how I can do that?


texts: Join based on maximum overlap in PostGIS/PostGresQL?

texts:
I have two sets of polygons in two tables. The sets overlap each other. For each polygon in set A, I would like to get the ID of the polygon in set B that it overlaps the most. I'm using PostgreSQL with the PostGIS extension.

I know just enough about SQL to know that you can only join based on true/false conditions. So this won't work:

SELECT
  a.id as a_id,
  b.id as b_id,
FROM
  a
JOIN
  b
ON
  max(ST_Area(ST_Intersection(a.geom, b.geom)))

because max() can't be in the ON clause.

ST_Intersects() is a true/false test, so I could join on that, but polygons in set A will often overlap with more than one polygon in set B, and I need to know which one overlaps the most. ST_Intersects would presumably just return the first overlapping ID it came across, regardless of the extent of the overlap.

This seems like it should be do-able, but it's beyond me. Any thoughts?


texts:
I'm trying to classify my study area by Semi-Automatic Classification Plugin

When I want to compare the spectral differences between two types of land-use, it shows that FloatingPointError: invalid value encountered in greater_equal

How could solve this error?


texts: Extracting polygon attributes from covering polygon with highest area using ArcPy

texts:
I have a fishnet grid (polygon layer) and another polygon layer containing some hexagons.  I would like for each square polygon in the fishnet layer that covers a hexagon (even a tiny edge of it) to take the name attribute of the hexagon.  If more than one hexagon touches a square polygon, I want the square polygon to take the name of the hexagon which covers the highest area within that square.

I’m thinking something like:

gp.SpatialJoin_analysis(Fishnet, Hexagons, Output1, "JOIN_ONE_TO_ONE", "KEEP_COMMON", "Name 'Name' true true false 12 Text 0 0 ,First,#,Hexagons.shp,Name,-1,-1", "INTERSECTS", "0 DecimalDegrees", "")

How can I replace “INTERSECTS” with a different rule i.e. hexagon with largest area?

A spatial join analysis may not be the best analysis here and I’m open to suggestions.  I’m using ArcGIS Desktop 9.3 and Python 2.5.


texts: Fundamental concepts of ESRI Javascript API

texts: I am building a Web Application for the first time using the ESRI's Javascript API and wanted to know the basic fundamental concepts that I should keep in mind?


texts: Checking whether a GeoServer service is WMS or WMTS?

texts:
I'm using GeoServer 2.10.1 connected to a PostGIS database. I have created a new workspace and activated the WMTS checkbox in the workspace definition. The service URL always contains wms either when I activate WMS or WMTS option. Then, I guess that this is not a good way to check whether the service is WMS or WMTS...

The service URL is as follow:

http://my.server.fr/geoserver/myworkspace/wms?

A layer complete URL is a follow (copied from Layer Preview page):

http://my.server.fr/geoserver/myworkspace/wms?service=WMS&version=1.1.0&request=GetMap&layers=test&styles=&bbox=1466568.5,6127476.88999999,1528548.733,6176365.98899999&width=768&height=605&srs=EPSG:3947&format=application/openlayers

I can't see any other place than the worspace definition where to specify this option for a given service (neither within stores or layers definition).

I need to make sure that the service is WMTS and not WMS. I have checked layers display within QGIS but I'm not able to find this information out within QGIS GUI.

Is there a simple way to check this out?


texts: FeatureToPolygon (to cut polygons by line) not working correctly outside of ArcGIS

texts:
I asked a question yesterday, about splitting polygons using a line, and I ended up using the approach that is described in the comments here.

Although this worked with my test data, it causes a lot of issues when running it with more complex datasets.

Therefore I am working on a new approach, which works fine in ArcGIS, but not outside of ArcGIS.

In the following screenshot you can see my polygons and my line. I would like to have each polygon, crossed by a line, to be split by that line (in this case, each polygon would turn into two polygons).

Here is the workflow:

Step 1) MakeFeatureLayer_management: Create feature layer (because we are selecting the polygons)
Step 2) SelectLayerByLocation_management: Select all the buffers that are crossed by the line, because only those buffers should be considered
Step 3) buffersToBeSplit: iterate over the selection to extract all names, and push them into a list (each buffer has a unique name, which is found in the column Name)
Step 4) SelectLayerByAttribute_management: iterate over the list, and use each name to select the corresponding buffer
Step 5) FeatureToPolygon_management: for each buffer: use the buffer to cut it with the line

Note: I am splitting each feature separately, because splitting all of the buffers at once creates a lot of extra polygons (see comments mentioned above)

And here is the code:

import arcpy

polygon = r"D:\Data\buffers.shp"
line = r"D:\Data\line.shp"

# Step 1
polygon_fl = arcpy.MakeFeatureLayer_management(polygon,'polygon_fl')

# Step 2
arcpy.SelectLayerByLocation_management(polygon_fl,'INTERSECT',line)

# Step 3
buffersToBeSplit = []
for row in arcpy.SearchCursor('polygon_fl'):
    buffersToBeSplit.append(str(row.Name))

# Step 4, 5
i = 0
for bufferName in buffersToBeSplit:
    i += 1
    query = '"Name" = \'%s\'' % bufferName
    arcpy.SelectLayerByAttribute_management('polygon_fl','NEW_SELECTION',query)
    arcpy.FeatureToPolygon_management(['polygon_fl',line],r"D:\Data\Output\cutPoly%s.shp" % str(i))

Proof that the Select by Location works:

And here is the result when I am running this code from the Python console in ArcGIS:

Exactly what I need! Each buffer has now been split into two seperate parts.

Now the problem arises when I am running this code from outside of ArcGIS (PyScripter). The result looks as follows:

You can see that only the very first polygon is split (top left) and all the others are ignored. The fact that each polygon is returned means that the the selection is working fine, but what is not working is the part that uses the line to cut that selection.

I've been trying for quite a while now, and I absolutely cannot figure out why this code would work inside of ArcGIS but not outside.

PS: what I am trying to implement is what can be easily done in SAGA GIS using the Polygon-line-intersection tool . If you are a QGIS-user, you can find the tool in the Processing toolbox.


texts:
I'm using OL3 with the searching API Nominatim. I have a search bar and I execute this :

  http://nominatim.openstreetmap.org/search/"+string+"?    format=json&addressdetails=1&limit=1&polygon_XXX=1

The API give multiple POLYGON format :

polygon_geojson=1 : Output geometry of results in geojson format.
polygon_kml=1 : Output geometry of results in kml format.
polygon_svg=1 : Output geometry of results in svg format.
polygon_text=1 : Output geometry of results as a WKT.

My aim : Draw a polygon, the city borders.

My problem : I didn't find an example to use these types of result, can I have help please ? I think I have to create a new layer to add to the map but I don't know how to create the layer with the resource.


texts: Assurance that ST_Envelope works with geography as expected

texts:
Using PostGIS, I have a table of data with geographic points defined as last_known_coordinates GEOGRAPHY(POINT,4326). I want to find all of these within a geographic bounding box with given top left and bottom right geographic coordinates, and be sure that the calculations are geographic, not geometric, e.g. that the results account for the curvature of the earth.
The query I'm using contains the fragment: last_known_coordinates && ST_MakeEnvelope(?, ?, ?, ?, 4326)::geography
As you can tell I'm a bit new to this, and I know what I want, but not necessarily how to do it correctly. Other answers hint that this is correct, but I'd like to know for sure: Will this point data and SQL fragment take into account the curvature of the Earth in the bounding box?


texts: Show map using esri leaflet, with 4326 spatial reference

texts:
I have this MapServer

http://geoservices.big.go.id/arcgis/rest/services/RBI/Rupabumi/MapServer

It turns out it uses 4326 spatial reference.

I want to present it using leaflet, my code so far:

    var map = new L.Map('map', {
    crs: L.CRS.EPSG4326
  }).setView([-6.1333, 106.75], 0);

  L.esri.tiledMapLayer('http://geoservices.big.go.id/arcgis/rest/services/RBI/Rupabumi/MapServer', {
    maxZoom: 14,
    minZoom: 0,
    continuousWorld: false,
    noWrap : false
  }).addTo(map);

But the map is not showing up. I found that http://geoservices.big.go.id/arcgis/rest/services/RBI/Rupabumi/MapServer/tile/0/20/25 will give me a tile, but the request from leaflet stops at /0/0/0. Can anybody help me?


texts:
I want to prepare the heat maps so what parameters are required for that??

I have large point data (I want to represent the density on the map)

I gone through the following posts :

How to build effective heat-maps?
http://forums.arcgis.com/threads/43069-How-to-create-a-heat-map-from-shapefile-data

I am using ArcGIS Desktop 10.1 and ArcGIS Server 10.1 + ArcGIS Java Script API 3.0


texts:
So I am a new user with GRASS 7.2, in Windows, and the problem was when I first booted GRASS GIS, I had to create a new location, for my data, but every time I was trying to create that location GRASS would either

Crash

or popup an empty error

Is there a way to resolve that, or do I need to install an older version of GRASS?
.


texts: Why can't QGIS 3.16 & KML Tools export point symbology for KML/KMZ

texts: I am trying to export styled points from QGIS to KML/KMZ format. The points are styled according to an attribute value ("Subsidence"). I want to export the KML/KMZ with the symbology, but when I open the result in Google Earth, it isn't preserved. Can someone help me understand why and if there is a way to accomplish this?


texts: Point shapefile not showing up when adding layer in QGIS?

texts:
The layer has all the points in the attribute table and when I click "Zoom to Layer" it zooms to the location where the points should be, but no points are appearing.

I also have a tried uploading the points as a CSV-file which did successfully create the points. Unfortunately, they showed up in Canada and they should be in Colorado.

I checked all the coordinate systems and made sure everything was lining up as it should and still nothing.


texts: QGIS: How to create a crosstab using SAGA toolbox cross calssification and tabulation?

texts: I am trying to use Cross classification and tabulation SAGA toolbox. I have a raster DEM file (14 classes) and raster file precipitation (13 classes). I want to do cross-tab - but I don't know how many classes I've to choose (look at the picture below)in SAGA toolbox. And when table is created I don't know which data layer is in rows and which in columns.


texts: Replacing cloud values with NA in a raster stack

texts:
I am removing clouds from a raster stack of 768 pixel_qa bands which are provided by USGS as part of surface reflectance product. The details of raster stack (qa) are:

> qa
class      : RasterStack
dimensions : 2520, 2714, 6839280, 768  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 205755, 287175, 3950955, 4026555  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Before I use this qa raster stack to mask NDVI stack, I use the following code to replace cloud values to NA:

> qa[qa %in% c(224, 480, 992)] <- NA

where 224 is cloud pixel value for Landsat 4-7 and 480 & 992 for Landsat 8.

The code runs forever! What could be wrong here?


texts: Iterate over two Feature Classes Inputs in ModelBuilder

texts:
I use a tool (Centrality of Urban Network Analysis) that requires two input files: a polygon layer made of buildings and a Network Dataset. Since my area of analysis is quite huge I divided it into 40 parts. I put all the polygon layers in a folder and all the ND into another.

I want to iterate the process with this tool but problem is, in ModelBuilder, we can only use one iterator, Iterate Features Classes in my case. Since we cannot use two iterators in ModelBuilder, I wonder how it is possible to iterate a process when someone is using a tool that requires two different inputs.


texts: Failing to connect to an authenticated WFS using owslib.wfs.WebFeatureService

texts:
Does anyone know why the following code isn't working? Substituting the exact same code works for the WMS at the same URL, so not sure what pyproj has to do with it (see traceback below).
I can access an authenticated WFS in QGIS, but connecting through owslib fails. I can access an unrelated and unauthenticated WFS using this code.
I intentionally haven't included the username and password for this post.
from owslib.wfs import WebFeatureService
bom_wfs = WebFeatureService(url="https://spatial.bom.gov.au/cgi-bin/mapserver/users/username/wxs?", username="username", password="password", version="1.0.0")
print(bom_wfs.identification.title)

The traceback error code I get is:
 C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\python.exe C:/Users/jmckechn/PycharmProjects/TreeProject/wfs_test.py
Traceback (most recent call last):
  File "C:/Users/jmckechn/PycharmProjects/TreeProject/wfs_test.py", line 8, in <module>
    bom_wfs = WebFeatureService(url="https://spatial.bom.gov.au/cgi-bin/mapserver/users/username/wxs?", username="username", password="password", version="1.0.0")
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\owslib\wfs.py", line 50, in WebFeatureService
    timeout=timeout, headers=headers, auth=auth)
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\owslib\feature\wfs100.py", line 104, in __new__
    auth=auth,
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\owslib\feature\wfs100.py", line 144, in __init__
    self._buildMetadata(parse_remote_metadata)
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\owslib\feature\wfs100.py", line 171, in _buildMetadata
    feature, featuretypelist, parse_remote_metadata, auth=self.auth
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\owslib\feature\wfs100.py", line 398, in __init__
    wgs84 = pyproj.Proj("epsg:4326")
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\proj.py", line 295, in __init__
    self.crs = CRS.from_user_input(projparams if projparams is not None else kwargs)
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\crs.py", line 226, in from_user_input
    return cls.from_string(value)
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\crs.py", line 189, in from_string
    return cls.from_epsg(val)
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\crs.py", line 167, in from_epsg
    return cls("epsg:{}".format(code))
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\crs.py", line 146, in __init__
    super(CRS, self).__init__(projstring)
  File "pyproj/_crs.pyx", line 307, in pyproj._crs._CRS.__init__
  File "C:\Users\jmckechn\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\lib\site-packages\pyproj\datadir.py", line 62, in get_data_dir
    "Valid PROJ.4 data directory not found."
pyproj.exceptions.DataDirError: Valid PROJ.4 data directory not found.Either set the path using the environmental variable PROJ_LIB or with `pyproj.datadir.set_data_dir`.

Process finished with exit code 1


texts: 1 year ArcGIS Desktop advanced student not working after formatting PC?

texts: I had ArcGIS 10.x installed in my PC and it was a student version for learning. I had to format the PC once and its just a month after installing ArcGIS 10.x but I cannot authorize it with the same code in the same machine. It but says its already authorized multiple times.


texts: arcpy.ImportMetadata_conversion works for file geodatabase but not SDE

texts:
I want to update metadata (the abstract/description) on an Esri ArcGIS feature class. I am using the Esri arcpy python lib. When I try this using a file geodatabase it works fine. When I try it using a SDE definition (SQL Server connection) it fails with the following error message:

  ExecuteError: Error using XSLT transformation "c:\program files
  (x86)\ArcGIS\desktop10.3\Metadata/Stylesheets/gpTools/merge imported
  metadata with existing.xslt": An error occurred while loading document
  '__TestData\myfeatureclass_xslttransforma3.xml'. See InnerException
  for a complete description of the error. Failed to execute (Merge
  Import With Existing). Failed to execute (ImportMetadata).

It fails on the line that imports the metadata into the layer. Here is the complete code:

import traceback
import os
import sys
import arcpy
import xml.etree.ElementTree as ET

def XSLTFilePath():
    installDir = arcpy.GetInstallInfo("desktop")["InstallDir"]
    return os.path.join(installDir, r"Metadata\Stylesheets\gpTools\exact copy of.xslt")

def metadataScratchFilePath(featureClassName):
    metadata_temp_folder = arcpy.env.scratchFolder
    return os.path.join(metadata_temp_folder, featureClassName + ".xml")

def nodeIsDefined(node):
    if node == None:
        return False
    elif len(node) == 0:
        return False
    return True

def updateMetaDataFileAbstract(metadata_file,abstractText):
    # update or add abstract info to the metadata
    doc = ET.parse(metadata_file)
    root_node = doc.getroot()
    dataIdInfo_node = root_node.find("dataIdInfo")
    if nodeIsDefined(dataIdInfo_node):
        abstract_node = dataIdInfo_node.find("idAbs")
    else:
        dataIdInfo_node = ET.SubElement(root_node, "dataIdInfo")
        abstract_node = ET.SubElement(dataIdInfo_node, "idAbs")
    abstract_node.text = abstractText
    tree = ET.ElementTree(root_node)
    tree.write(metadata_file)

def updateFeatureClassAbstract(featureClassLocation,featureClassName,abstractText):
    # set up scratch file location for XML
    metadata_file = metadataScratchFilePath(featureClassName)
    if os.path.exists(metadata_file):
        os.remove(metadata_file)
    # create XML scratch file
    arcpy.env.workspace = featureClassLocation
    arcpy.XSLTransform_conversion(featureClassName, XSLTFilePath(), metadata_file)
    # update XML file with new abstract
    updateMetaDataFileAbstract(metadata_file,abstractText)
    # write updates to the feature class
    arcpy.ImportMetadata_conversion(metadata_file, "FROM_ARCGIS", featureClassName,
                                    Enable_automatic_updates=False)

try:
    print("This works...")
    USERPROFILE = os.getenv('USERPROFILE')
    WS_FILE = USERPROFILE + '\\Documents\\ArcGIS\\Default.gdb'
    updateFeatureClassAbstract(WS_FILE,"myfeatureclass","This is my new abstract text")
    print("This doesn't...")
    updateFeatureClassAbstract('My_DB.sde', "myfeatureclass", "This is my new abstract text")
    # normal completion
except:
    crash_info = traceback.format_exc()
    crash_msg = "Script Crashed."
    print(crash_msg)
    print(crash_info)
    print(crash_msg)
    print(crash_info)
    sys.exit(99)
    # abnormal completion