texts: Polyline changes into polygon in google map drawing tool

texts:
I have been using google map provided tools where it can draw the polygon and polyline accordingly. I will first click on the polyline icon and start plotting and suddenly when I double click it turn to be polygon and I notice icon which I clicked earlier have been not in selected mode now. Below is my codes where I do the calling when I do the overlay completed. So where could be my error ?

google.maps.event.addListener(drawingManager, 'overlaycomplete', function(e) {
    drawingManager.setDrawingMode(null);
     var newShape = e.overlay;
            newShape.type = e.type;
     if(e.type.toString()=="polygon"){
        //customised info window
     }
     else if(e.type.toString()=="polyline")
     {
       //customised info window
     }
     else if(e.type.toString()=="marker")
     {
      //customised info window
     }
 });
google.maps.event.addListener(drawingManager, 'drawingmode_changed', clearSelection);


texts: How to delete overlay within the same feature class?

texts:
So, the following scenario:
I have plenty of polygons that are within another polygon, all in the same feature class. All of the polygons (the one overlaping the other ones as well as the ones being overlapped) have attributes which are important for my project. What I basically want is to erase those parts of the selected polygon (see screenshot below) that overlap with other polygons in the same feature class.

http://www.bilder-upload.eu/upload/8440e2-1461147192.jpg

This is basically just reparing an error which went unnoticed for some time during an editing session and using the auto-complete-polygon option. The selected feature is supposed to be railtracks... you can surely guess what happened there. First step would be just to delete the overlaps, the rest I can do manually with the split-polygon tool.

//Edit: I'm using Esri ArcGIS 10.2.2


texts: Plotting elevation maps and shaded relief images from latitude, longitude, and elevation pairs

texts:
I have a set of latitude, longitude, and elevation pairs (roughly a grid of such values, though it is not uniform), and I'd like to be able to plot an elevation map and perhaps also a shaded relief image for this data. And I'd like to do this using python.

While I've seen some interesting examples of people making such maps online, they always start by reading in some data file which is in some format which I am not familiar with (perhaps they have a measurement device which simply outputs this format).

I'm curious how I might go about doing such a thing from the simple latitude, longitude, and elevation pairs.

I currently only need to plot small geographic regions (10 square miles, for example).


texts: How to populate auxiliary fields in QGIS 3?

texts:
After creating a user field in the auxiliary storage, how can I populate this field with a default value? Normally I'd do this in the attribute table with "Update All" or "Update Selected", but the auxiliary fields aren't available in the pulldown menu. The same is true in Field Calculator.

I've also tried Select All, and then Edit->Modify Attributes of Selected Features. I can enter the default value in this dialogue, but when I click OK the fields do not actually change. I get the same result if I try to enter a default value in the Attribute Form editor in Layer Properties (and I've selected "Apply default value on update").


texts: MapInfo Pro 15.2 incorrect coordinates when updating column

texts: I have an issue getting the correct coordinates when updating the columns, if I double click on the individual object then it comes up with the coordinates in UTM. However, even when I change the map projection to UTM then the coordinates remain the same. I have tried to save the files as a copy and attempt again, but unsure what to do now.


texts:
I am trying to clip a series of lines selecting only the geometry portion outside of a polygon using PostGIS.

From my reading on what st_intersection does i.e. Returns a geometry that represents the shared portion of geomA and geomB. I'd like to do the opposite of this so return the geometry that isn't shared by geomA and geomB.

I've been trying to work with st_difference which should return a geometry that represents that part of geometry A that does not intersect with geometry B. But i'm having some trouble getting the result I want. This is what i've written so far...

SELECT
    t.id, ST_Difference(t.the_geom_webmercator, t1.dock) as intersection
FROM
    trips t,
    (select
        d.id, d.name,ST_Buffer(ST_Centroid(the_geom_webmercator),1000) as dock
    from
        docks d) as t1
WHERE
    ST_Intersects(t.the_geom_webmercator,t1.dock)
    t.id IN (25558, 25555, 25553, 25533, 25517, 25523, 25501, 25506,25532,25439)
    and ST_Difference(t.the_geom_webmercator, t1.dock) is not null

This returns

the pink lines are the result of the above query, so some lines stop at the circle boundaries and others don't
the circles are the polygons I want to clip everything by,
the underlying point line is the result of basically the same query
as above replacing st_intersection with st_difference which
returns the inverse of the desired result i.e. the geometries within
the circles.

Any suggestions on how to amend this query?


texts: Make a leaflet TileLayer appear above a GeoJSON layer

texts:
I'm trying to add a layer of place labels above a GeoJSON layer, so that these labels remain visible. However when I add the layer per the code below, the tileLayer remains below the GeoJSON layer.

    var CartoDB_PositronOnlyLabels = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_only_labels/{z}/{x}/{y}.png', {
        attribution: '&copy; OpenStreetMap &copy; CartoDB',
        subdomains: 'abcd',
        maxZoom: 19,
        zIndex: 3
    });
    CartoDB_PositronOnlyLabels.addTo(map);

Adding a few more things that don't work for the future reader:

According to this leaflet PR you could add a pane: 'overlayPane' option upon layer creation.
The hack appearing here and here only works if the tileLayer you're adding is of tileJSON format

var topPane = map._createPane('leaflet-top-pane', map.getPanes().mapPane);
var topLayer = L.mapbox.tileLayer('lxbarth.map-vtt23b1i').addTo(map);
topPane.appendChild(topLayer.getContainer());
topLayer.setZIndex(9);


texts: How do I get a non-circle marker in a category map on CartoDB?

texts:
I have a CartoDB map which has two layers. I want both to display, and I want both to be category maps — there are 3 categories in each which display as three different colors. But I want to distinguish between them by marker shape. Ideally, one layer will be a circle and the other will be a triangle or a square.

I see that I can choose "ellipse" or "arrow" as a marker-type in CartoCSS. Arrow is not useful for my purposes as it doesn't have enough space to show color. Apparently "rectangle" was supposed to exist, but hasn't been implemented for non-torque maps.

Chris from CartoDB suggested using an SVG icon, which DOES allow being filled in with the given marker-fill. I'm guessing this is because the fill for the square in the SVG file is "Undefined," and thus able to be inherited. Unfortunately, when I do this, the squares look terrible on my map because they're all pixelated. They also don't have sharp borders, even though I've provided a border in the image.

The only solution I've been able to think of thus far is to create PNG/SVG icons with the required colors, upload them and use them directly for each category. This is obviously not dynamic, but at least it would be aesthetically acceptable.

Is there some other solution I'm not thinking of? Does anyone know why the SVG seems to be rendered and then resized such that it becomes pixelated?


texts:
I have a polygon shapefile and a (merged) line shapefile (both are zipped).  When I use Vector -> Analysis Tools -> Sum Line Lengths on these shapefiles, I get a layer with some of the LENGTHS filled with numeric values. However, a significant portion of the LENGTHS are now 'nan', even when the COUNT is a non-zero value.

Here are a few possible explanations, but I'm not sure which one is correct:

Error with the shapefiles.
Shapefiles are too big.
The data type of the LENGTH field is inappropriate.
I need to unzip them.

I've tried doing the same with a few of the unmerged line shapefiles and this problem still seems to occur.


texts: interpolate raster into a continuous surface image--QGIS

texts:
I have a raster of non-continuous surface with valid values and nodata values, and I'm trying to interpolate it into a continuous surface image without nodata values using QGIS. Anyone could help?
The non-continuous surface (white area has nodata value and coloured area has valid values:

What I want to achieve using QGIS:

Thank you!


texts: How to join data from NHGIS to make choropleth maps in R/ggplot2?

texts:
I have downloaded shapefiles of county/state boundaries and census data keyed to those shapefiles from the NHGIS. Both the shapefiles and the census data have a GISJOIN column to be used for merging the data. I'm able to do this just fine in QGIS, but I can't for the life of me get it to work in R/ggplot2.

This is what the data looks like:

str(churches_1850)
'data.frame':   36 obs. of  30 variables:
 $ GISJOIN : Factor w/ 36 levels "G010","G050",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ YEAR    : int  1850 1850 1850 1850 1850 1850 1850 1850 1850 1850 ...
 $ STATE   : Factor w/ 36 levels "Alabama","Arkansas",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ STATEA  : int  10 50 60 90 100 110 120 130 170 180 ...
 $ COUNTY  : logi  NA NA NA NA NA NA ...
 $ COUNTYA : logi  NA NA NA NA NA NA ...
 $ AREANAME: Factor w/ 36 levels "Alabama","Arkansas",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ AET001  : int  579 114 1 114 12 6 56 885 282 430 ...
 $ AET002  : int  18 NA NA 4 NA NA NA 5 69 187 ...
 $ AET003  : int  NA NA NA 252 NA NA NA 1 46 2 ...
 $ AET004  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ AET005  : int  NA NA NA NA NA NA NA NA 2 5 ...
 $ AET006  : int  17 2 1 101 21 8 10 20 27 24 ...
 $ AET007  : int  5 1 NA 1 NA NA 1 6 2 10 ...
 $ AET008  : int  NA NA NA 5 9 1 NA 2 6 89 ...
 $ AET009  : int  NA NA NA NA NA NA NA NA 3 5 ...
 $ AET010  : int  NA NA NA 2 NA NA NA NA NA NA ...
 $ AET011  : int  1 NA NA NA NA 2 NA 8 42 63 ...
 $ AET012  : int  NA NA NA NA NA NA NA NA NA NA ...
 $ AET013  : int  577 168 5 185 106 16 87 809 405 779 ...
 $ AET014  : int  NA NA NA NA NA NA NA 1 2 57 ...
 $ AET015  : int  162 52 3 17 26 6 16 97 206 282 ...
 $ AET016  : int  NA NA NA NA NA NA 2 NA NA NA ...
 $ AET017  : int  5 7 18 12 3 6 5 8 59 63 ...
 $ AET018  : int  NA NA NA NA NA NA NA NA 2 NA ...
 $ AET019  : int  NA NA NA NA NA NA NA NA 4 5 ...
 $ AET020  : int  4 5 NA 4 1 NA NA 16 30 5 ...
 $ AET021  : int  1 NA NA 5 NA 1 NA NA 4 1 ...
 $ AET022  : int  3 NA NA 22 NA NA NA 3 7 15 ...
 $ AET023  : int  3 13 NA 10 2 NA NA 7 25 13 ...

This is what I'm trying so far:

library(rgdal)
library(ggplot2)
states_1850 <-readOGR(dsn = "shp/nhgis0007_shapefile_tl2000_us_state_1850/",
                        layer = "US_state_1850")
churches_1850 <- read.csv("csv/nhgis0007_ds10_1850_state.csv")
states_1850_df <- fortify(states_1850, region = "GISJOIN")
states_with_data <- merge(states_1850_df, churches_1850, by.x="id", by.y="GISJOIN")

ggplot() +
  geom_polygon(data=states_with_data,
               aes(x=long, y=lat, group=group, fill=AET001)) +
  coord_map()

I've tried a number of other ways of merging the data, but they all either end up not merging correctly, or tearing the shapes of the states/counties, or being so large that plotting them crashes my machine (8GB ram, shouldn't be a problem).

This seems like it should be easier. What's the proper way to do this?

UPDATE: This method of merging data frames seems to be the older way of doing things. Using geom_map is much more powerful, since it merges the data on the fly rather than at the beginning. It also allows for easier faceting. See the ggplot2 docs for geom_map.

My revised plotting code looks like this:

p <- ggplot() +
  geom_map(data = churches_1850,
           aes(map_id = GISJOIN, fill = AET015),
           map = shp_1860, size=0.01, alpha=.9) +
  expand_limits(x = shp_1860$long, y=shp_1860$lat) +
  scale_fill_gradient2(low="white", high="red")
print(p)


texts:
I am looking to create a Python script that will take X,Y coordinates from one field in an attribute table, compare them to another set of X,Y coordinates in a different field, and populate a field with the shortest distance between those two points based on a road layer.  So just to visualize it..

Origin:X,Y      Destination:X,Y         Distance: 130.34km


texts: Projecting gps data to get length of resulting line shape in ArcMap 10.0

texts:
I have gps data in WGS84 geographic coordinates.  I will use the gps data to build a line shape in ArcMap 10.0.  The data falls in UTM 15N.

To get most accurate length of line, should I project the point data to the UTM and then build the line in the UTM.  Or should I build the line in WGS84 then project the line data to the UTM.  Or does it not matter?  What are the effects?

If UTM data doesn't best preserve distance then I will add that UTM must be used.


texts: Creating custom legend to map to custom css torque heatmap

texts:
I'm hoping someone can help me with html to create a legend that maps to a custom heatmap I made using the style CSS. When I try using the default legend, it neither maps to the right custom colors, nor does it have values from my torque-aggregate-function.

Picture of map + legend here: https://cl.ly/1l1v062Q3K11

Custom CSS here:

Map {
  -torque-frame-count: 1;
  -torque-animation-duration: 30;
  -torque-time-attribute: "cartodb_id";
  -torque-aggregation-function: "sum(total_visits)";
  -torque-resolution: 3;
  -torque-data-aggregation: linear;
}
#layer {
  marker-width: 12;
  marker-fill: white;
  marker-fill-opacity: 1;
  marker-file: url(https://cartodb-libs.global.ssl.fastly.net/cartodbui/assets/unversioned/images/alphamarker.png);
  marker-allow-overlap: true;
  marker-line-width: 1;
  marker-line-color: #FFFFFF;
  marker-line-opacity: 1;
  image-filters: colorize-alpha(blue, cyan, lightgreen, yellow , orange, red);
}


texts: draw multipart polygon with one surrounding border

texts:
As you can see on the image below, the waterbody is divided by the border line because it consists of several polygons.

I could merge the polygon, but I dont want to change the shapefile.
Is there a option in QGIS to draw te waterbody just with a border sourrounding its "absolute extend"?


texts: Different results for same transformation on different PostgreSQL/Postgis installations


texts: How to have a date earlier than 1753 in QGIS date field?

texts:
I'm trying to map out some of my family tree but I've stumbled upon a problem of not being able to put in early dates.

I created the shapefile in QGIS and included a date field (dd-mm-yyyy).
Whenever I enter in a date earlier than 1753 It will not let me (i.e. it wont let me put a 6 next to the 1 for 1600) it then defaults back to todays date.

How do I get around this??


texts: Running ArcGIS API for JavaScript locally?

texts:
I'm pretty new to arcgis javascript api, and I'm wondering how to run it locally.
I have downloaded the api, and put the folder "arcgis_js_api" in the root folder of my apachi server, and edited the files "init.js" and "dojo.js" and replace the strings '[HOSTNAME_AND_PATH_TO_JSAPI]dojo' to 'localhost/arcgis_js_api/library/3.12/3.12/' as said in the documentation. But still it won't work. am I missing something?


texts: I am trying to develop an accuracy assessment tool using ArcObjects. As an output I would like to show a matrix with the fields such as road ,land,water etc.  I just want to know how can I fetch unique values from the  raster value attribute table into a matrix format.


texts:
I'm looking for an algorithm to calculate the 3d bounding box of a spherical polygon
 - i.e. the min/max X, Y and Z over that part of the surface of a unit sphere.

The purpose is to aid fast area-overlap calculations for large numbers of simple (convex) spherical polygons :
Having all the bounding boxes will let you quickly calculate, for a given target polygon, a large majority of the other polygons which definitely do not intersect it (because their bounding boxes are disjoint).

I already have a calculation that does this type of selection, using polygon centres and maximum-radii.
However, that will take O(N^2) time to calculate (for all source polygons; for all target polygons).
Whereas, this approach would hopefully run in O(N.log N) time, as the XYZ information can be ordered and then sorted/searched (e.g. using kdtrees).

I think this is a bit less trivial than it sounds.
For instance, the max X value on a polygon boundary segment (which is a great circle) can obviously occur at a point somewhere between the two endpoints.
But also, in fact, an extreme value can occur somewhere inside the polygon away from the boundary
-- for example, we must be able to get a max Z=1.0 for a polygon drawn around the North pole, even though its boundary lies entirely South of the pole.

Can anyone help with this ?


texts: Polyline geometry to graphic - SimpleLineSymbol is the only type not displaying on the map

texts:
I have modified the Measurement widget to keep a session based history of the measurements the user has made. When a user clicks on a history item, it should display the geometry associated with that history item as a GraphicsLayer on the map. I am using knockout to manage the history items and to retrieve measurement metadata when a history item is clicked.

At this point, both my Polygons (for area) and Points (for location) work just fine with the SimpleFillSymbol() and the SimpleMarkerSymbol(), respectively. However, the Polyline geometry returned from a distance measurement is not displaying on the map with the SimpleLineSymbol().

Here's the code:

var graphicLayerId = "measurementHistoryGraphicsLayer";

function addGraphicsLayerToMap(graphicsLayer) {
    var lay = getGraphicsLayerFromMap();
    if (lay !== undefined) {
        lay.clear();
        lay.add(graphicsLayer);
        map.removeLayer(lay);
    }
    map.addLayer(graphicsLayer);
}

function createGraphicFromGeometry(viewModel) {
    //Determine the symbol type
    var symbol;
    switch (viewModel.activeTool) {
        case "area":
            symbol = new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID,
                                                      new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new esri.Color([255, 0, 0]), 2),
                                                      new esri.Color([255, 0, 0, 0.25]));
            break;
        case "distance":
            symbol = new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new esri.Color([255, 0, 0]), 3);
            break;
        case "location":
            symbol = new esri.symbol.SimpleMarkerSymbol(esri.symbol.SimpleMarkerSymbol.STYLE_SQUARE, 10,
                                                            new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new esri.Color([255, 0, 0]), 3),
                                                            new esri.Color([0, 255, 0, 0.25]));
            break;
    }
    var graphic = new esri.Graphic(viewModel.geometry, symbol, { "extent": viewModel.extent, "unitName": viewModel.unitName });
    return graphic;
}

function createGraphicsLayerFromGraphic(graphic) {
    var graphicLayer = new esri.layers.GraphicsLayer({ id: graphicLayerId });
    graphicLayer.add(graphic);
    graphicLayer.setRenderer(new esri.renderer.SimpleRenderer(graphic.symbol));
    return graphicLayer;
}

function getGraphicsLayerFromMap() {
    return map.getLayer(graphicLayerId);
}

$(document).on('click', '#emv_measurement_history .list-group-item', function () {
    $('#emv_measurement_history .list-group-item.list-group-item-info').removeClass('list-group-item-info');
    $(this).addClass('list-group-item-info');
    var measurementData = ko.mapping.toJS(ko.dataFor($(this)[0]));
    var graphic = createGraphicFromGeometry(measurementData);
    var graphicsLayer = createGraphicsLayerFromGraphic(graphic);
    addGraphicsLayerToMap(graphicsLayer);
    map.setExtent(measurementData.extent);
});

Like I said, this works fine for both area and location, but distance does not seem to work. I've even tried adding a hard-coded polyline value in there are creating a SimpleLineSymbol from that without success.

For additional information, here is the Polyline info:

[
    [
        [
            2591769.2297164765,
            5236836.417134136
        ],
        [
            2573584.2281166334,
            4620357.96034264
        ],
        [
            2557384.1428811993,
            4038303.8136230526
        ],
        [
            3124973.8484519687,
            4260007.60486125
        ],
        [
            3714518.451309448,
            4454862.77067183
        ],
        [
            4324318.833989203,
            4618552.510359674
        ],
        [
            4666465.839330839,
            4693607.843734423
        ],
        [
            5013294.285789721,
            4757423.375729576
        ]
    ]
]

And the spatial reference is set to 102100.


texts: Python Shapely could not find 'geos_c.dll' in Windows 10

texts:
I was trying to install and work on 'Shapely' library of python on Windows 10. Installation was successful but when I tried to import it then it gave error :

OSError: Could not find lib geos_c.dll or load any of its variants []

Steps Followed :
pip3 install shapely

(Successfully installed shapely-1.7.0)
Then I imported library in python3 shell :
from shapely.geometry import Point

And got the error :

OSError: Could not find lib geos_c.dll or load any of its variants [].

Platform:

Windows 10
Python 3.8.1


texts: Installing PostGIS on older homebrew versions of PostgreSQL (e.g. postgresql@9.6)?

texts:
As there are two ways to solve this question, I asked for another approach here.

Whenever a major release of PostgreSQL appears, homebrew will happily install it right away, unless one pins a certain version or installs a specific version from Homebrew/core.

So far so good, but the problem is that PostGIS also keeps updating and eventually will be incompatible with pinned or Homebrew/core versions of PostgreSQL. Again, one can pin PostGIS, yet Homebrew/core does not offer previous versions. From time to time I do brew cleanup and somehow I am now left with a broken PostGIS database. Ok, so "time to jump to PostgreSQL 10.1", I though, but pg_upgrade now gets stuck with the error

  ERROR: could not access file "$libdir/postgis-2.3": No such file or directory

That means I can't back nor forth: I cannot figure out how to get PostGIS working with postgresql@9.6 and I cannot pg_upgrade to the new version of PostgreSQL (10.1).

Therefore my question can be answered in two ways. Here I am asking:

How to install PostGIS with older PostgreSQL versions, i.e. postgresql@9.6, using homebrew?

The second approach is described at Upgrading PostgreSQL with broken PostGIS installation ("$libdir/postgis-2.3" error)?

Tried so far:

This gist offers a solution and it might actually work. However, I had to use ln -s instead of cp -a. It seems to be working now, but unfortunately I have the problem that my database is already on 2.4 and there seems to be no downgrade path. Will keep trying.


texts: Command line join attributes by location using OGR

texts:
I have two overlapping vector files (polygons and sample) and want to do a spatial join from the command line to extract all features in polygons that intersect with sample. Essentially I want to have the same result as if I was using the "Join attributes by location" function in QGIS with "Only keep matching records" including the attribute table.

This should be possible with OGR and a SQL query, but I can't get the syntax to work.
So far I have tried with a vrt file following this description

And then something like:

ogrinfo -sql "SELECT * FROM sample WHERE ST_INTERSECTS (polygons,sample)" -dialect SQLITE plogonesandsample.vrt

However I find this very cumbersome and not generic.

I seek a solution that works from the command line.


texts: PostGIS: ERROR: lwpoly_from_lwlines: shell must have at least 4 points

texts:
I have achieved to close lines in order to build a polygon with:

st_addpoint(
    st_makeline(geom ORDER BY .cd_noseq),
    st_pointn(
        st_makeline(geom ORDER BY cd_noseq),
    1)
)

but then, when I try to make a polygon with st_makepolygon like this:

st_makepolygon(
    st_addpoint(
        st_makeline(geom ORDER BY cd_noseq),
        st_pointn(
            st_makeline(geom ORDER BY cd_noseq),
        1)
    )
) as geom

I get the following error:

ERROR:  lwpoly_from_lwlines: shell must have at least 4 points
********** Error **********

ERROR: lwpoly_from_lwlines: shell must have at least 4 points
SQL state: XX000

Technology:

  PostgreSQL 9.2.7 on x86_64-unknown-linux-gnu, compiled by gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-4), 64-bit"
  "POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.6" LIBJSON="UNKNOWN" (core procs from "2.1.0 r11822" need upgrade) TOPOLOGY (topology procs from "2.1.0 r11822" need upg (…)


texts: In an attribute table I am making new fields bounded each other by the Field Calculator. Once I have written the formula for the new field in the Field Calculator, ArcGIS make the calculation for the all field. If, in the field calculator I say that Field B = Field A + 1, and I calculate Field B, then I change value to Field A, I notice that Field B does not consequently change values. I would need to type the formula for Field B again to refresh its value. How can I make field B dependent on Field A? I have also used Ilwis software and there was just a box to be checked in its Field Calculator window; is there in ArcGIS something similar?


texts: Topology processing after osm2pgrouting?

texts:
I imported an OSM dataset using osm2pgrouting, and checked the topology like this:

pgr_sfbay=# select pgr_analyzegraph('ways', 0.001, id:='gid');
NOTICE:  PROCESSING:
NOTICE:  pgr_analyzeGraph('ways',0.001,'the_geom','gid','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...
NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 2045
NOTICE:                            Dead ends: 115909
NOTICE:  Potential gaps found near dead ends: 99024
NOTICE:               Intersections detected: 11525
NOTICE:                      Ring geometries: 15
 pgr_analyzegraph
------------------
 OK
(1 row)

Then I tried pgr_nodeNetwork() and pgr_createTopology() on the imported table 'ways' and found that the number of isolated segments and gaps are significantly reduced:

pgr_sfbay=# select pgr_analyzegraph('ways_noded', 0.001);
NOTICE:  PROCESSING:
NOTICE:  pgr_analyzeGraph('ways_noded',0.001,'the_geom','id','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...
NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 320
NOTICE:                            Dead ends: 27912
NOTICE:  Potential gaps found near dead ends: 11343
NOTICE:               Intersections detected: 47398
NOTICE:                      Ring geometries: 15
 pgr_analyzegraph
------------------
 OK
(1 row)

So I have a few questions:

Why osm2pgrouting gives different results than pgrouting? Does it use the same functions? Or is it because it has some info like overpass that are in the original OSM data?
Is it a good idea to run pgrouting topology functions after osm2pgrouting importing? If so, how to handle the cost, reverse_cost, one_way columns, etc?

For this test I am using mapconfig.xml from GitHub without any changes.


texts: Creating GetCapabilities document in more languages

texts:
I'm providing INSPIRE metadata and services and I want them in more languages than one. I expect that most users are locals, but some users are definitely from other European countries. Therefore I want to support national language and English. In metadata documents, I have used gmd:PT_FreeText_PropertyType, but GetCapabilities file follows other XSD schemas (WMS or WFS, according to the service).

Do you know about some possibilities of adding additional language to GetCapabilities document?


texts: ArcGIS Field Calculator passes European-formatted floats as tuples

texts:
I have the problem that the ArcGIS Field Calculator doesn't read a float number when it is included in Python by referring to the field that contains it, in this way:

def giveResult(floatValue):
  if floatValue>0:
    a=(1000,5000,10000,15000,20000)
    for i in range(len(a)):
      if a[i]/1000>floatValue:
        return a[i]
        break
      else:
        return 3
  else:
    return 0

This code returns only the values 3 and 0, because the line

if a[i]/1000>floatValue:

doesn't work properly. The problem seems to be my European Windows localization settings, where decimals are displayed with commas instead of points. These settings are applied to float fields in ArcGIS. It seems that Field Calculator interprets these floats as tuples because of the comma.

How can this data be handled?

I finally have it, in this way: checking which is the variable type in each case, and only if the type is tuple, we catch the first part (it's enough to build the function); we don't mind if the variable isn't tuple (don't need to catch that condition).

Only that part of the code would be:

  for i in range(len(a)):
      if type(OFEhu)==tuple:
        OFEhu=OFEhu[0]
      if (a[i]/1000)>OFEhu:
        return a[i]


texts: Determine the feature dataset of a layer in an Esri file geodatabase

texts:
I have an Esri file geodatabase with a single feature dataset.
Using ogr and python
ds = ogr.Open('myfile.gdb')

for i in range(ds.GetLayerCount()):
    layer = ds.GetLayerByIndex(i)

I can iterate the layers of the GDB file.
While using ogr and python how can I determine which feature dataset the layer belongs to?

Looking at the GDAL/OGR API I haven't found anything.
The driver description only mentions how to create a layer in a specific feature dataset.
This has an idea, but I hope for a more simple way to do it.


texts: how to write Topology Graph Edges geometry to polyline feature class

texts:
i have created topologygraph of a polygon feature . and iwant individual ids of polylines

i want to write Topology Graph Edges geometry to new polyline feature class

any samples OR documentation about this task?


texts: Getting bad allocation while exporting using QGIS?

texts:
I am using QGIS 3.4.

I'm creating a pretty simple map that shows major roads, cities/towns, and state/county boundaries for 4 northeastern U.S. states. All but the boundaries have labels (not buffered). I was successfully exporting large PDFs (22x36, 300 dpi) of this map as I progressed to make sure my file wasn't too large. At a certain point, my map stopped exporting, so I started turning layers off to see what would and wouldn't export. Now I'm at a point where I'm lucky if I can export just the basemap at 75 dpi. My maps renders at the same speed, but I get the 'bad allocation' notice within seconds of attempting to export anything.

Any suggestions?


texts: arcpy CalculateFieldManagement - single quotes in expression

texts:
I'm using the following to write a value to a field.

arcpy.CalculateField_management (newPha, "PHACODE", """'{}'""".format(name), "PYTHON_9.3")

The only problem is that the string I'm attempting to write to the table has an apostrophe. (ex. "Saint Anne's Bay").

The single quote is clearly throwing off the expression, but the triple-double quote solution isn't working.
Any ideas on how to deal with this???


texts: Can I create areas where all label are forced away

texts: I have some points overlaying a raster. These points have some quite complicated labels which block out the raster so I would like them to sit outside the edge of the raster and point at the point to which they refer. Is it possible (with or without maplex) to set the labels no not touch an area such as a raster or a feature in another layer? (I have the outline of the raster as a polygon)


texts: Permission Issues on ArcGIS Server with SOC machines

texts:
I am trying to figure out some permission issues for ArcGIS Server. I have 2 SOC machines that are associated with SOM. I noticed this problem when I added these SOC machines to SOM.
This service that I have published is a GP service.

Appreciate any ideas.

Thanks
Jay

First Error:
~~~~~~~~
WARNING Can not open directory c:\arcgisserver\arcgisoutput. The system cannot find the path specified. Probable cause: The network path is inaccessible to ArcGIS Server Object Container account. 2011-11-08 11:14:15 SDM.ServerDirectoryManager cwf414

Second Error:
~~~~~~~~~~~
ERROR Can not open file D:\Projects\GPServices.mxd. The system cannot find the path specified. Probable cause: The file is inaccessible to Server. 2011-11-08 11:07:31 PRJ/GPServices.MapServer cwf415

ERROR Server Object instance creation failed on machine cwf415. Can not open file D:\Projects\GPServices.mxd. The system cannot find the path specified. Probable cause: The file is inaccessible to Server.  2011-11-08 11:07:31 PRJ/GPServices.MapServer cwf415

ERROR Server Object instance creation failed on machine cwf415. Can not open file D:\Projects\GPServices.mxd. The system cannot find the path specified. Probable cause: The file is inaccessible to Server.  2011-11-08 11:07:31 PRJ/GPServices.MapServer cwf415


texts: Installing 1-band color raster color table?

texts:
I would like to use the "1-band raster colour table" extension.
I'm using QGIS Wroclaw 1.7.4.

In the plugin installer I can't find this extension, even I quote for third party extensions thing. Does anybody know its name in the Plugin manager? Is there another plugin similar to 1-band raster that I could use?

I've read on another forum that this extension don't work on QGIS 1.7.... is that right?

bccaddRcolor: seems to be a similar extension than 1-band raster, but i work only on QGIS 1.8.

Where can I download it ?


texts:
I'm looking at putting together a light Python install on a thumb drive that I can take to client offices and do some basic GIS processing. I'm thinking of installing the following:

PyScripter
Python 2.6/2.7
Numpy
Scipy
GDAL/OGR with the Python bindings
GRASS, likewise with the Pythong bindings
QGIS

Does anyone have any suggestions of any libraries / software that might be a useful addition? Typically I'll be working with raster datasets, but any and all suggestions would be appreciated.

(I'll cw the question if people think that its appropriate).


texts: How can I get the timeslider using Cartojs and Torquejs

texts: I want to build the timesilder on my webpage.The timeslider on Carto Builder is nice for me.Now I can get the first img' timeslider, But I want the second one and don't know how to achieve it.What library or ui should I need?


texts: Scratch files left behind by arcpy.sa.Con

texts:
I have a Con() statement that creates two "extra" scratch files during it's processing, but that do not get cleaned up once the operation is over and the output is saved.

    r=arcpy.sa.Con(elv, rCtypeAlpine, rCtype, 'VALUE > 10000')
    r.save('covtypeprms_a')

prior to saving r, I can see that there is also a scratch file for that object. This gets renamed with r.save(). But I'm still left w/two directories (grids?).

I've been finding and deleting them using arcpy.ListRasters('t_t*'), which works, but this is really just a hack. Lock files seem to get left behind with this approach. I've also tried del r to get rid of the handle, but not destroy the copy on disk ('covtypeprms_a').

Why are these appearing and is there more precise way to reference these intermediate objects (so's I can delete them!)?

update on enviroment

product: Desktop 10.1 sp1
OS: Win7 Enterprise sp1

python: Enthought 32-bit in use via Wing IDE

sys.executable yields 'C:\Users\joeblow\AppData\Local\Enthought\Canopy32\User\Scripts\python.exe'


texts: Range value of raster data varies before and after mosaic implementation

texts:
I am trying to raster::mosaic bunch of .tif files extracted from hdf4 MODIS files. My script is working fine but I am puzzled by the range of raster values (min, max) before and after data processing. The original hdf4 rasters are 8 bits integer raster ranged between 0-255 for every individual raster files before mosaicing. However, when I do a raster mosaic then the resulting raster is ranged between 0-25. My question is why?

I tried to reproduce this through a reproducible example but the problem doesn't apply to any example data such as below:

#------------------------
rx1 <- raster(nrows=2400, ncols=2400, vals = floor(runif(5760000, min=0, max=256)) )
extent(rx1) <- c(11119505, 12231456, -4447802, -3335852)
rx2 <-  raster(nrows=2400, ncols=2400, vals = floor(runif(5760000, min=0, max=256)) )
extent(rx2) <- c(12231456, 13343406, -4447802, -3335852)

lr.1 <- list()
lr.1[[1]] <- rx1
lr.1[[2]] <- rx2

lr.1$fun <- mean
rast.mosaic1 <- do.call(mosaic,lr.1)
rast.mosaic1

#class       : RasterLayer
#dimensions  : 2400, 4800, 11520000  (nrow, ncol, ncell)
#resolution  : 463.3129, 463.3125  (x, y)
#extent      : 11119505, 13343407, -4447802, -3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names       : layer
#values      : 0, 255  (min, max)

If the same line of codes is applied to this data, then the result is:

lr.2 <- list()
lr.2[[1]] <- r[[1]]
lr.2[[2]] <- r[[3]]

lr.2$fun <- mean
rast.mosaic2 <- do.call(mosaic,lr.2)
rast.mosaic2

#class       : RasterLayer
#dimensions  : 2400, 4800, 11520000  (nrow, ncol, ncell)
#resolution  : 463.3127, 463.3127  (x, y)
#extent      : 11119505, 13343406, -4447802, -3335852  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
#data source : in memory
#names       : layer
#values      : 0, 25.4  (min, max)


texts: Performing multiple selection using MapBasic?

texts:
I am new to MapBasic programming.

I am trying to make multiple selection from a table where the ID number exist in my array and save all the values into another table which called results.
ARRContainsValue() is a function which checks the values. This function works well. My only problem is the table result only contains the last selection, not all of them. Can anyone help me with this.

fetch first from structurepoint
for i=1 to tableinfo(structurepoint,8)
    if ARRContainsValue(Route_IDs_Used(),route.id)=1 then
        select * From structurepoint where rowid=i into results noselect
    end if
    fetch next from structurepoint
next


texts: How to handle time-based raster data in the ArcGIS Server JS API

texts:
The TimeSlider functionality in the ArcGIS Server JS API is great for working with time-based vector features - you can scroll forwards and backwards through time to see the vector features which existed at a given time instant.

Has anyone created a similar time-slider for rasters?

I have a series of monthly weather grids which I'd like to display, allowing the user to scroll to any month (eg, show me the average rainfall for June).

The time-slider relies on the featureclass being defined in ArcMap with a Time attribute, but this isn't available for raster layers.

The best workaround I can think of is:

create a separate map document for each raster
build a separate tiled map service from each map document
add a DOJO slider to the map, with increments to mark each raster
use JavaScript to switch the appropriate raster on at each mark on the DOJO slider

Is anyone aware of a better/faster approach? This will be quite a laborious process as it will mean creating and maintaining multiple MXDs and MSDs.

Thanks


texts: Reclassifying floating point raster in ArcGIS Desktop?

texts:
I have a floating point raster having dimensions of 478 x 418 columns and rows respectively. The cell size is 30 x 30 meters, therefore, it covers an area of 179.5 Sqkm.

The intent is to obtain one row for each pixel of the map in a table.

Is there any way around except 'Raster to Points' tool available in ArcMap?

As such, I would like to retain values as pixels to compute other statistics required for validation of this raster through Success Rate and Prediction Rate Curves! Like, frequency of training/testing points in each cell.

How should I reclassify it by retaining its original values in contrast to assigning new values using 'Reclassify tool?


texts:
I've created a basic map using D3 with countries GeoJSON. Here is the demo.

Now when the user clicks on any coordinate on the map, I am showing a weather info within a tooltip with weather icon as a marker.

countries = countriesGroup
    .selectAll("path")
    .data(json.features)
    .enter()
    .append("path")
    .attr("d", path)
    .attr("id", function(d, i) {
        return "country" + d.properties.iso_a3;
    })
    .attr("class", "country")
    // add a mouseover action to show name label for feature/country
    .on("mouseover", function(d, i) {
        d3.select("#countryLabel" + d.properties.iso_a3).style("display", "block");
    })
    .on("mouseout", function(d, i) {
        d3.select("#countryLabel" + d.properties.iso_a3).style("display", "none");
    })
    // add an onclick action to zoom into clicked country
    .on("click", function(d, i) {
        var eventLocation = d3.mouse(this);
        var coordinates = projection.invert(eventLocation);
        var proxy = "https://cors-anywhere.herokuapp.com/";
        var wetherInfoUrl =
            "https://api.darksky.net/forecast/c68e9aaf0d467528b9363e383bde6254/" + coordinates[1] + "," + coordinates[0] + "?exclude=minutely,hourly,daily";
        $.ajax({
            url: proxy + wetherInfoUrl,
            success: function(response) {
                tooltipDiv
                    .transition()
                    .duration(200)
                    .style("opacity", 0.9);
                tooltipDiv
                    .html('
Dynamic Weather info: '+ response.currently.humidity +'

')
                    .style("left", (eventLocation[0] - 250) + "px")
                    .style("top", (eventLocation[1] - 100) + "px");
            }
        });
        d3.selectAll(".country").classed("country-on", false);
        d3.select(this).classed("country-on", true);
        boxZoom(path.bounds(d), path.centroid(d), 20);
    });

Now the issue is, the tooltip remains at the absolute position relative to the body (It stays at the exact position on the screen which is misleading if we pan and zoom the map) whereas I want the tooltip to be relative to the coordinate clicked on the map and should be fixed to the clicked coordinate no matter where we zoom or pan in the map, similar to this Leaflet example (In short I want tooltip to behave similarly to the typical map pin markers).


texts: imposm: Failing to import OSM data to PostgreSQL

texts:
Somehow I can't import OSM data to my database. I'm using Debian, PostgreSQL 9.1 and PostGIS 1.5. The user root has been granted all privileges for the database osm.

Input:

imposm -U root -d osm -m imposm-mapping.py --read --write --optimize sweden-latest.osm.pbf

Output:

Enabling Shapely speedups.
loading imposm-mapping.py as mapping
password for root at localhost:
[16:35:06] ## reading sweden-latest.osm.pbf
[16:35:07] coords: 11724k nodes: 234k ways: 1674k relations: 11k (estimated)
[16:39:56] coords: 17332k nodes: 75k ways: 1445k relations: 7k
[16:39:57] reading took 4 m 51s
[16:39:57] ## dropping/creating tables
Traceback (most recent call last):
  File "/usr/bin/imposm", line 9, in <module>
    load_entry_point('imposm==2.4.0', 'console_scripts', 'imposm')()
  File "/usr/lib/python2.7/dist-packages/imposm/app.py", line 231, in main
    db.create_tables(tag_mapping.mappings)
  File "/usr/lib/python2.7/dist-packages/imposm/db/postgis.py", line 129, in create_tables
    self.create_table(mapping)
  File "/usr/lib/python2.7/dist-packages/imposm/db/postgis.py", line 160, in create_table
    pg_geometry_type=mapping.geom_type))
psycopg2.ProgrammingError: permission denied for relation spatial_ref_sys
CONTEXT:  SQL statement "SELECT SRID         FROM spatial_ref_sys WHERE SRID = new_srid"
PL/pgSQL function "addgeometrycolumn" line 75 at SQL statement
SQL statement "SELECT AddGeometryColumn('',$1,$2,$3,$4,$5,$6)"
PL/pgSQL function "addgeometrycolumn" line 5 at SQL statement

EDIT:

And this is what I get when I use the user postgres instead:

[13:36:35] ## dropping/creating tables
[13:36:36] ## writing data
Traceback (most recent call last):
  File "/usr/bin/imposm", line 9, in <module>
    load_entry_point('imposm==2.4.0', 'console_scripts', 'imposm')()
  File "/usr/lib/python2.7/dist-packages/imposm/app.py", line 243, in main
    writer.relations()
  File "/usr/lib/python2.7/dist-packages/imposm/writer.py", line 52, in relation                                    s
    cache = self.cache.relations_cache()
  File "/usr/lib/python2.7/dist-packages/imposm/cache/osm.py", line 59, in relat                                    ions_cache
    return self._x_cache(self.relations_fname, RelationDB, mode, estimated_recor                                    ds)
  File "/usr/lib/python2.7/dist-packages/imposm/cache/osm.py", line 68, in _x_ca                                    che
    cache = x_class(x, mode, estimated_records=estimated_records)
  File "tc.pyx", line 104, in imposm.cache.tc.BDB.__init__ (imposm/cache/tc.c:12                                    63)
IOError: 3


texts: CPUs to process lidar data using ArcGIS for Desktop?

texts:
I am looking to build a computer to process data in ArcGIS 10.2.2.

The main use will be to rasterize 10-100's GBs of lidar data.

However, I am unclear as to the capabilities of ArcGIS 10.2.2 to utilize hyper-threading and multiple cores. In prior iterations of ArcGIS, I remember that only some geoprocessing tools could utilize this added computational capability. I am therefore trying to decide between the Intel Core i5-4690K and i7-4690K.

Does anyone have any experince with these processors and ArcGIS 10.2.2? Even better, does anyone have experience with these cores for lidar data processing?


texts: How to assign polygon attributes to grid cells using centroids in QGIS?

texts:
I am relatively new using QGIS.

Unfortunately I can´t figure out the following yet.
I have two polygon shapefiles, one is a grid I have created using QGIS 2.12.3-Lyon (Creat graticule) and the second represents native communities´ lands. I want to assign to each graticule´s cell the name and attributes of the corresponding overlapping community. Some cells overlap with more than one communal land. In these cases I want to assign the name of the community that overlaps with the cell centroid.

How can I do this?


texts: R MODIS MOD13Q1 - How to handle data quality issues?

texts:
I'm running into some issues with 250m 16-day MOD13Q1 NDVI data using the MODIS package in R.

Looking at the following raster plot, there are certain days where the NDVI does not match the true surface characteristics. Specifically ndvi.2017.11.01 and ndvi.2018.01.17. Even ndvi.2018.03.06 and ndvi.2018.02.18 are questionable. Is it possible to programmatically exclude these layers? Is there something to gain by using the following sds layer: MODIS_Grid_16DAY_250m_500m_VI:250m 16 days VI Quality?

Below is the code to produce the above plot:

library(raster)
library(MODIS)
library(parallel)
library(pbapply)
library(magrittr)
library(dplyr)
library(stringr)

# Final Study Area Boundary
study_area <- readRDS(gzcon(url('http://web.pdx.edu/~porlando/study_area.RDS')))
# Study CRS
epsg_26910 <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "

# Download Data
runGdal(product = "MOD13Q1", collection = "006"
           ,tileH = 9, tileV = 4
           ,begin = str_replace_all("2017-09-01", "-", ".")
           ,end = str_replace_all("2018-09-30", "-", ".")
           ,overwrite = TRUE
           ,job="temporalComposite"
           #,SDSstring="" # not sure about which SDS string to use...
)

# Process data
ndvi_path <- "./MODIS/MOD13Q1.006" # 250 m

ndvi_files <- list.files(path = ndvi_path, pattern = "\\.hdf$"
                     ,all.files = FALSE, full.names = TRUE
                     ,recursive = TRUE
                     ,ignore.case = FALSE)

processNDVI <- function(file_path, study_area = study_area, proj = epsg_26910) {

  cat("\n")

  # extract date from file path
  date <- stringr::str_extract(file_path, '[0-9][0-9][0-9][0-9].[0-9][0-9].[0-9][0-9]')
  date_clean <- stringr::str_replace_all(date, "\\.", "-")

  # coerce 16 day average to daily composite?
  dates <- seq.Date(from = as.Date(date_clean)
                    #,to = as.Date(date_clean) + lubridate::days(15)
                    ,to = as.Date(date_clean)
                    ,by = "1 day")

  sds <- get_subdatasets(file_path)

  ndvi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days NDVI", sds)] %>% readGDAL %>% raster
  evi <- sds[grepl("16DAY_250m_500m_VI:250m 16 days EVI", sds)] %>% readGDAL %>% raster
  # VI Quality may be useful?
  vi_quality <- sds[grepl("16DAY_250m_500m_VI:250m 16 days VI Quality", sds)] %>% readGDAL %>% raster

  # combine NDVI and EVI into a single stack
  r <- stack(ndvi, evi)
  names(r) <- c(paste0("ndvi.", date), paste0("evi.", date))

  r <- raster::projectRaster(from = r, res = 250, method = 'ngb', crs = proj)
  study_area <- spTransform(study_area, CRSobj = proj)

  if(identical(crs(study_area), crs(r))) {

    m <- mask(r, study_area)
    cr <- crop(m, study_area)

    pblapply(dates, function(x) {
      pblapply(1:nlayers(cr), function(y) {
        layer_name <- names(cr[[y]])
        var <- gsub("\\..*$", "", layer_name)
        date <- gsub("-", "\\.", x)
        writeRaster(cr[[y]]
                    ,filename = paste0("./output/", var, "/", var, ".", date, ".tif")
                    ,format = "GTiff"
                    ,overwrite = TRUE)
      })
    })
}

pblapply(ndvi_files, function(x) {
  processNDVI(file_path = x
              ,study_area = study_area
              )
  }
  ,cl = detectCores()-1
    )

ndvi <- stack(list.files(path = "./output/ndvi/", pattern = ".tif$", full.names = T))
plot(ndvi)


texts: Why esriSpatialRelTouches relation does not work?

texts:
I have two feature layers. The first one is just a polygon and the other one consists of a series of polygons which are enclosed by the first one. Thus, some polygons of the second feature layer touch on the boundaries of the first feauture layer. I wrotre a code (the basic part shown below to identify which polygons of the second layer touches the boundary of the polygon in the first layer but while the code runs without a problem it does not identify (e.g. by giving FID with Msgbox) any such polygon whilsty it should.
Any help please?

Dim pBlockBoundaryFC  As IFeatureClass
Set pBlockBoundaryFC = pBlockBoundary.FeatureClass

Dim pBlockBoundaryFields  As IFields
Set pBlockBoundaryFields = pBlockBoundaryFC.Fields

Dim pGeom As IGeometry

Dim pSpatialFilter As ISpatialFilter
Set pSpatialFilter = New SpatialFilter

Dim pBlockBoundaryCursor As IFeatureCursor
Set pBlockBoundaryCursor = pBlockBoundaryFC.Search(pSpatialFilter, False)

Set pBlockBoundaryFeature = pBlockBoundaryCursor.NextFeature
Do Until pBlockBoundaryFeature Is Nothing

Set pGeom = pBlockBoundaryFeature.Shape
With pSpatialFilter

'**************
'Conversion
Dim pPolygonBoundary As ITopologicalOperator
Set pPolygonBoundary = pFeature.Shape

Set pGeom = pPolygonBoundary.Boundary
'*************

Set .Geometry = pGeom
    .GeometryField = "Shape"
    .SpatialRel = esriSpatialRelTouches
End With

Dim pRelOp As IRelationalOperator
Set pRelOp = pGeom

  Set pNewParcelsFeature = pNewParcelsCursor.NextFeature
  Do Until pNewParcelsFeature Is Nothing

    FID = pNewParcelsFeature.Value(intposFID)
    If pRelOp.Touches(pNewParcelsFeature.Shape) Then

    ''Report the FID of each parcel that 'touches' the exetrnal boundary
    MsgBox FID
    End If

Set pNewParcelsFeature = pNewParcelsCursor.NextFeature

Loop

Set pBlockBoundaryFeature = pBlockBoundaryCursor.NextFeature

Loop


texts: How to calculate the area of overlap of a buffer in ArcMap

texts: I am a beginner with GIS. I have used the buffer tool in ArcMap to draw a 1 km buffer around a set of three survey posts. I want to calculate the area of where the buffer's overlap. Is there a tool in ArcMap which will allow me to do this?


texts: Georeferencing an oblique angle photograph in repeat-photograph studies

texts:
Backround: Aerial photography and remote sensing can be used to measure current land cover and recent land use change. However, these technologies are fairly recent and therefore time series are limited when using them. In contrast, historical landscape photos exist for many locations from much further back in time. Repeat-photography studies return to the locations of historical photographs a take identical photographs to examine changes. However, in general, these studies are limited to qualitative measures of change inferred manually because these obliquely taken photographs can't easily be georeferenced.

Question: I'm helping with an repeat-photography study in Alaska, we have historical (up to 100 years old) photos of the landscape (with lat/long) and we'll be retaking photos this year at the same locations. We also have access to a high resolution DEM and aerial imagery. The goal is to quantify forest cover change between the photos. I'm looking for a method to georeference these obliquely taken photographs, drape them over the landscape, and create a map to be compared with aerial photography. Essentially I want to transform an oblique photo to an orthophoto.


texts: Cannot select input layers for the "Polygonize" tool in QGIS

texts:
I am trying to convert a raster layer to a polygon layer using QGIS 3.16. 2. When I open the polygonize tool, the program does not allow me to select any existing raster layers from the 'Input layer' drop-down menu.
I have tried manually selecting the raster from my computer (a .TIF tile), but just get the following message:

The Log Messages window does't seem to offer any further information, at least in the OGR tab, reading:
2021-01-16T11:42:49     WARNING    Cannot open C:/Users/User/Desktop/GBG_slope.tif.()

I am not sure what's happening here. What potential ways can I solve this issue?


texts: Using the processing framework in a PyQGIS standalone application

texts:
I asked a question several days ago about creating an in-memory layer using the processing.runalg() command in a standalone PyQGIS application (Performing in-memory difference operation between vector layers in PyQGIS).  The solution seems to work perfectly well in the QGIS Python console, but not when I try it in a standalone application.

The relevant snippet of code I am running is:

diffp = processing.runalg("qgis:fixeddistancebuffer", layer, 0.000001, 10, False, None)
diff = processing.getObject(diffp['OUTPUT'])

(layer is a vector layer defining line segments)

My desired output is this, and this is what I get when I run the commands from the Python console in QGIS:

But what I get from the vector layer I create in the standalone application is shown in green in the following image (I've zoomed in so that it's visible):

What I get from the standalone application appears to be only one feature from my desired output.  I'm running the same two lines of code for both the console and standalone applications, and there are no exceptions thrown by the code.  How can I produce the entire desired result from a standalone application, rather than just a single feature?

EDIT:

I've tried using the following line of code to directly create a .shp file, in case the in-memory layer was somehow causing the problem:

diffp = general.runalg("qgis:fixeddistancebuffer", layer, 0.000001, 10, False, "shape_file.shp")

This produces the same results as when I use the in-memory layer, so the problem I'm experiencing is definitely coming from runalg().


texts: QGIS modeler - how to extract polygons by multiple arguments?

texts: I'm using QGIS modeler to extract lines by not one, but multiple attributes. For example, I have a road network, and I want to extract only road with character values in field TYPE: value_one and value_two. How to do this using QGIS modeler?


texts: How to locate a house on the 3D Map correctly?

texts:
I've done my 3D map generation in WebGL/THREE.JS project and have an issue with the building locations.

I have the data about building position ( its polygon ) from DB server, but first of all I want just to locate not the polygon data, but a simple cube. The aim is clear is as a crystal: to set a building to the correct position by putting a cube on the 3D map.

Now it does look like:

I have demonstrated it, because I want to admit, that to create a building and place it on the map - isn't a problem. The problem is to set building into the correct position/location on the map.

What I have done before?

Because the building record has two fields in DB:

position on map ( lat | long ) | polygon data of building ( points collection )

I thought that I must convert the geospatial coordinates of each building: latitude, longitude to its cartesian values, because WebGL/THREE.JS don't work with the shperical/geographic coordinates system.

So I have built also a JSON web-service which helps in data transformations and similar aims.

For e.g. the conversion of the point { 38.23, 55.89 } to its cartesian representation:

Now... I'm really confused how to use it for locating the building correctly on the map.
If you would like see, I'm prinring the position of each building in console:

Spherical coordinates and converting at web-service side by such formula in my C# code:

internal Models.PointCartesian ConvertFromSphericalToCartesian(double angleTheta, double anglePhi)
{
    var radTheta = ConvertDegreesToRadians(angleTheta);
    var radPhi = ConvertDegreesToRadians(anglePhi);

    var model = new Models.PointCartesian()
    {
        X = earthRadius * Math.Cos(radTheta) * Math.Cos(radPhi),
        Y = earthRadius * Math.Cos(radTheta) * Math.Sin(radPhi),
        Z = earthRadius * Math.Sin(radTheta)
    };

    return model;
}

So... Right now I'm thinking about some proportions to use... But this way can be incorrect... So what's your suggestion how to set a cube into the correct position on the 3D map just for showing the real position of some buildings.

PS:
The records lat | long of each building in DB is 100% correct.
Thanks!


texts: Run kriging in ModelBuilder with various shapefiles

texts:
I need to run kriging in a model, using 5 different point shapefiles. I put as input %Name% to run the different shapefiles, but in Z value field I have an error, because it's different in each shapefile.

So, how do I run kriging through various shapefiles that have different Z value fields? I want 5 output rasters. I'm doing this:


texts: Arcpy: MakeFeatureLayer_management: Where is the file?

texts:
I am trying to understand arcpy/python, and am confused about where a layer file is located during the process.

This is my code:

import os, sys
import arcpy

gdb  = "C:\\Example.gdb"
featureClass = "C:\\Example.gdb\feature"
arcpy.MakeFeatureLayer_management(featureClass, "OutLayer", "\"Field_Name\" = 'Name One'")
arpy.CopyFeatures_management("Outlayer", "C:\\Example.gdb\\test")

This seems to work, but I don't understand where "Outlayer" exists. If I run the code again, I get the message that I can't because "Outlayer" already exists. It is not in the .gdb, or the default.gdb. When I specify a path (e.g. C:\GIS\OutLayer) it isn't there. Finally, when I search for it using Windows Explorer, I can't find the file. It is not important, but just trying to figure out what is going on.


texts: How to use Geokettle to load shapefiles into oracle spatial

texts:
I'm looking for an example how to load a shapefile from the file system into an oracle spatial database.

What I have so far is:

I load the shapefile with the 'Shapefile File input' step
I use the calculator step to transform the WKT from the shapefile into a Geometry
I use the Table output to write the data into the database.

I get an error from the Table output step which says 'Error setting value #6 [GEOM Geometry] on prepared statement (Geometry), java.sql.SQLException: invalid datatype

The target column is of type SDO_GEOMETRY.

Thanks
D3


texts: Downloading OSM data for multiple polygons at one time

texts:
I am interested in the number of buildings found within each protected area in Africa. I have a protected area layer (WGS84) that I am using. The amount of data per protected area will be very small, most likely less than 200 nodes per protected area.

However, I cannot find a GIS tool (either QGIS or ArcGIS) that will allow me to download OSM data for multiple polygons from a vector layer at once, despite the data being very small. I did see that ArcGIS has an option under the OSM editor “Extent of data download” dropbox for intersection or union of inputs but these do not seem to work.


texts: QGIS: How to customize the list of file formats available when loading?

texts:
QGIS is able to load an impressive list of file formats, both vector and raster.  For example, here's a screenshot of just some of the raster formats available when using the Data Source Manager popup window at 3.10.5:

However, this eye-watering list is too extensive for my needs. I would like to remove those formats that I will never use, so that such popups contain only those formats that are of interest to me (for example, those highlighted in red).

I looked through Settings > Interface Customization, but did not see anything related.

My question: How to customize the list of file formats available when loading?


texts: Orthogonal projection of a point onto linestring in PostGIS

texts:
I am trying to find the orthogonal projection of a Point (in 4326 SRID) onto a LineString (in 4326). The picture shows the objective.

I have tried to use ST_Distance function, but in this case the function will only return the minimum distance to the starting point of LineString.


texts:
Does anyone know how to set a baseline elevation in Pix4D, or is that something that has to happen on the drone-capturing side?

The reason is that I have an area which is flown with LiDAR, several years have passed and now I want to measure the growth of the vegetation. However, I know the bare-earth baseline elevation for the take-off point, but the GPS vertical on the UAV is not that accurate, so if it is recording elevation with each picture that elevation may be relatively accurate for that flight, but could be ~5-10m off the baseline. If I could systematically adjust all elevations on the pictures by the same amount the drone is out from the bare earth DEM at the takeoff point then I would have a (hopefully) more accurate measurement of the difference between the LiDAR DSM and the UAV flown DSM.

Anyone know if that should happen when the UAV is set to fly, or if it is something that I can adjust in Pix4D before the processing?


texts: Converting coordinates from meters to decimal degrees in QGIS?

texts:
I have a shapefile. The coordinates (Lat/Long) are projected as meter (555000, 155000). But, I want to re-project it to decimal degrees(51.523666, -0.07467).

How can I do it with QGIS?


texts: Function execution hangs if connection to Internet is lost

texts:
I am trying to update feature layer. I am updating the records in a loop.

The problem is in the case Internet connection lost. It get stuck in the edit_features() function. I have to ctrl+c to stop script execution.

Is there a way to jump  to except: and continue the flow?

Here is my code:

from arcgis.gis import GIS
from IPython.display import display
gis = GIS("portal url", 'username', 'password')

#search for the feature layer named Ports along west coast
search_result = gis.content.search('title:Ports along west coast')

#access the item's feature layers
ports_item = search_result[0]
ports_layers = ports_item.layers

#query all the features and display it on a map
ports_fset = ports_layers[0].query() #an empty query string will return all
ports_flayer = ports_layers[0]
ports_features = ports_fset.features

# select San Francisco feature
sfo_feature = [f for f in ports_features if f.attributes['port_name']=='SAN FRANCISCO'][0]
sfo_feature.attributes
try:
    update_result = ports_flayer.edit_features(updates=[sfo_edit])
except:
    pass


texts: Vector data and geometry changes inexplicably when you zoom in

texts:
I am trying to make use of a multipolygon vector file representing electric utility service territories, with really interesting data tables. However when I read it in QGIS or OpenLayers it has a behavior I have never seen and can't control, which is that the data and geometry change as you zoon in. This is the file:
https://hifld-geoplatform.opendata.arcgis.com/datasets/c4fd0b01c2544a2f83440dab292f0980_0
And you can also examine the file as a WFS service at this address:
https://services1.arcgis.com/Hp6G80Pky0om7QvQ/arcgis/rest/services/Electric_Retail_Service_Territories/FeatureServer/0
From what I understand after much reading, the problem is the rendering order of the features. They should be rendered in descending order of size, so that the smallest is rendered last. In QGIS I can go Properties > Symbology > Layer Rendering > Control feature rendering order, and there I can enter "$area DESC NULLS LAST", which will change the way the file displays in QGIS. But that only changes the order in QGIS, whereas I want to upload this to PostGIS/Geoserver and display it with OpenLayers. So the question is:
In a multipolygon geojson, how can I reorder the features from largest to smallest?


texts: Stripping spatial reference from raster using ArcGIS Desktop?

texts:
I am creating an exercise for remote sensing students and would like to take an existing raster and remove the spatial reference from the file so that the students can practice using another image to orthorectify the first. I am using ArcGIS 10.

Any hints?


texts: Transport layer on OpenLayers 3 map: ol.layer.OSM is not a constructor

texts:
I'm using the Thunderforest Transport map and I want to overlay it onto my OpenLayers 3 map. After looking through the documentation and following their example (http://www.thunderforest.com/tutorials/openlayers/), I get the following error:

ol.layer.OSM is not a constructor

Everything in my code is working fine, but the following section raises this error:

var transport = new ol.layer.OSM("TransportMap", ["https://tile.thunderforest.com/transport/0/0/0.png?apikey=MY_API_KEY"]);
map.addLayer(transport);

I looked at the documentation (http://openlayers.org/en/latest/apidoc/ol.layer.html) and it seems that OSM is no longer present in ol.layer. What method do I use to add such a layer to my existing map? Should I add it as a Tile?


texts: Making NoData transparent in existing RGBA GeoTIFF

texts:
I have RGBA GeoTIFF where some areas have NoData values (255 255 255 255). Mapnik doesn't understand NoData and renders it as opaque white. How to make NoData areas as fully transparent?
I tried gdalwarp -dstalpha in.tif out.tif but it made NoData areas opaque black.


texts: PostGIS ST_EuclideanDistance() not working. Wrong upperleftx and y?

texts:
I'm trying to make this ST_EuclideanDistance() function work, but no success so far. I'm not getting the expected result because the raster seems to be drawn at a wrong location.

Documentation:
https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools

Firstly, I'm considering this test page to test the function:
https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/test

When running the second test case:

CREATE TABLE test_source_raster_1 AS
SELECT
    1 AS rid, rast
  FROM
    (SELECT
       ST_MakeEmptyRaster(10,10,0,0,1,1,0,0,4326) AS rast
     ) foo;
CREATE TABLE test_source_geometry_1 AS
SELECT
    generate_series(1,10) AS id,
    ST_RandomPoints(the_geom,10) AS the_geom
  FROM
    (SELECT
       ST_SetSRID(ST_Extent(rast::geometry), 4326) AS the_geom
    FROM test_source_raster_1) foo;
CREATE TABLE test_ref_raster_1 AS
SELECT
  1 AS rid, rast
  FROM
    (SELECT
       ST_MakeEmptyRaster(10,10,0,0,1,1,0,0,4326) AS rast
    ) foo;
 CREATE TABLE test_eudist_1_a_i_1 AS
   (SELECT
      1 AS rid,
      ST_EuclideanDistance(
         rast,'32BF','public',
         'test_source_geometry_1', 'the_geom', True
       ) AS rast
     FROM test_ref_raster_1);

I'm getting this result (every raster pixel is null):

And this is the expected result I should have gotten:


texts: pgr_dijkstra working but pgr_ksp is returning empty an resultset

texts:
This pgr_ksp function is not working, as answered in an earlier question

SELECT seq , id1 , id2 , cost FROM pgr_ksp(
    'SELECT gid as id, fnode_ as source, tnode_ as target, length as cost
     FROM road', 71670, 97745, 2, false) ;

I have tried all these solutions, but nothing works. Also there are no directions indicating possible restrictions that I could find.

pgr_dijkastra is working well with output.

SELECT seq, id1 AS node, id2 AS edge, cost,geom
FROM
   pgr_dijkstra('SELECT gid as id, fnode_ as source, tnode_ as target, length as cost
                 FROM road',71670 , 97745,false, false)
 JOIN road ON pgr_dijkstra.id2 = road.gid ;


texts: How to show QGIS plugin dialog always on top?

texts: Is it possible to define for a pyqgis plugin, that the plugin dialog is always shown on top, after starting it from the toolbar?


texts:
I made a web app that displays delivery routes associated with a user-input zip code. My question is, how can I get the map to zoom to the results? Here is my code. Im thinking of some sort of zoom-to-extent, is that along the correct path?

function init() {
myMap = new esri.Map("mapDiv", {
    basemap: "streets",
    center: [-89.9, 35.1],
    zoom: 10
});

dojo.connect(myMap, 'onLoad', function () {
    var scalebar = new esri.dijit.Scalebar({
        map: myMap,
        scalebarUnit: "english"
    });
});

require(['dojo/on', 'dojo/query'],
    function (on, query) {
        on(query('#myList'), 'change', function (e) {
            gpURL = e.target.value;
            console.log('Environment changed: ' + '"' + gpURL + '"');
        });
    });

if (!Array.prototype.indexOf) {
    Array.prototype.indexOf = function (obj, start) {
        for (var i = (start || 0), j = this.length; i < j; i++) {
            if (this[i] === obj) {
                return i;
            }
        }
        return -1;
    };
}

resultsLayer = new esri.layers.GraphicsLayer();
myMap.addLayer(resultsLayer)
}

dojo.addOnLoad(init);

function selectZIP(Zip_code) {

gp = new esri.tasks.Geoprocessor(gpURL);
gp.setOutputSpatialReference(myMap.spatialReference);

var params = {
    "ZIP": Zip_code,
    "Rte_Box": "C",
};

gp.execute(params, displayResults);

function displayResults(results, messages) {
    var featureset = results[0].value;

    var simplePolySymbol = new esri.symbol.SimpleFillSymbol();
    simplePolySymbol.setOutline(new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([0, 96, 170, .75]), 1));

    dojo.forEach(featureset.features, function (feature) {
        feature.setSymbol(simplePolySymbol);

        resultsLayer.add(feature);

    });

    myMap.addLayer(resultsLayer)
    myMap.Search.handleStats(results, messages)
}

myMap.Search = {
    gpGetStats: new esri.tasks.Geoprocessor(gpURL),
    getRouteStats: function () {
        console.log(dojo.byId('resultsLayer').value)
        var params = {
            "ZIP_CRID": dojo.byId('resultsLayer').value
        };
        resultsLayer.Search.gpGetStats.execute(params, resultsLayer.Search.handleStats);
    },

    handleStats: function (results, messages) {
        console.log(results);

        var content = '';

        content += '
';
        content += '
 ZIP CRID 	 Bus 	 Res 	 Tot 	 &lt;200 	 DS Key
';

        dojo.forEach(results[0].value.features, function (crid, i) {
            console.log(crid.attributes.ZIP_CRID + ', ');
            content += '
' + crid.attributes.ZIP_CRID;
            content += '	 ' + crid.attributes.BUS_CNT;
            content += '	 ' + crid.attributes.RES_CNT;
            content += '	 ' + crid.attributes.TOT_CNT;
            content += '	 ' + crid.attributes.LT_200_IND;
            content += '	 ' + crid.attributes.DS_KEY + '
';
        });
        content += '
';
        dojo.byId('details').innerHTML = content;
    },
}
}

dojo.ready(function () {
var clearMe = new dijit.form.Button({
    label: "Clear Map",
    onclick: function (init) {
        dojo.byId('details').innerHTML = '';
        resultsLayer.clear();
    }
}, "clearMe");
});


texts: JAXA RIKEN NowCast data reading in Python

texts:
I'm trying to read the RIKEN nowcast data in Python. They have a raw .dat file having a format of
4-byte float
plain binary,
little-endian

which contains 13 data layers ( coverage : from 60 N - 60 S ) with grid of [3600 x 1200 x 13]
import numpy as np
rnc_data = np.fromfile('gsmap_rnc_20200719.0700.dat', dtype = float)

print(rnc_data.min(), rnc_data.max(), rnc_data.shape)

The rnc_data array gives a shape of (28080000,) . What is confusing is (3600*1200*13)/28080000 gives exactly 2. Does this mean there's half the data there is supposed to be?
Am I missing something here?
Also if I had to convert this raw data to a raster, what would be the right way to do so


texts: Is there a non SDE version for the DisconnectUser (arcpy) function?

texts: My work environment consists of multiple file geodatabases, on a number of shared network drives that can be accessed by multiple users at a given time.  When I try to push out updates to feature classes or feature datasets, I get "locked by user" on many occasions.  I do not have SDE so was wondering if there is a way to accomplish a non SDE version of the DisconnectUser function so I can continue my updates.


texts: Filter a GeoPandas dataframe for points within a specific country

texts:
I have a dataframe of health survey data which is geotagged, called momdata. I have converted the dataframe to a GeoPandas dataframe, and used the geotags as the geometry column within the dataframe, for example:

momdata['geometry']

0        POINT (-0.0893 51.4735)
1       POINT (-0.0894 51.4732)
2         POINT (-0.0898 51.4717)
3        POINT (-0.0907 51.4727)
4       POINT (-0.0901 51.4723)
5       POINT (-0.0816 51.4742)

I want to filter the dataframe so that only points in the UK are returned. I have the UK coordinates in a separate Geoseries, called uk_geom, which I simply took from the GeoPandas built-in world map:

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
uk_index = world[world.name == "United Kingdom"].index
uk_geom = world.loc[uk_index, 'geometry']

This returns uk_geom as a GeoPandas GeoSeries:

type(uk_geom) # geopandas.geoseries.GeoSeries

I simply want to filter momdata for all points within momdata.geometry that fall within uk_geom, returning a dataframe of only UK-based survey observations. This should be simple.

I tried:

uk_momdata = momdata[momdata.geometry.within(uk_geom)]

but this returns an empty dataframe, when I know that some of the survey observations for sure are in the UK.

For example:

p1 = momdata.geometry.loc[0]
p2 = momdata.geometry.loc[1]

print(uk_geom.contains(p1)) # returns TRUE
print(uk_geom.contains(p2)) # returns TRUE

I tried this the other way round, checking which momdata points are contained within uk_geom:

uk_geom.contains(momdata.geometry).value_counts() #2040 false points

Also when I test the 'within' function on a point that I know is within momdata, I get an error:

print(point8.within(uk_geom)) # AttributeError: 'GeoSeries' object has no attribute '_geom'

I have assigned the correct co-ordinate reference system:

assert uk_geom.crs == momdata.crs # no problem

I also tried a basic 'apply' function using a predicate, but this returns an error:

momdata[momdata.geometry.apply(lambda p: uk_geom.contains(p))] # Null geometry supports no operations

I also tried a spatial join, but then I get the error that one of the join columns is not a DataFrame, as of course it's trying to join on the geometry column:

from geopandas.tools import sjoin
join_left_df = sjoin(momdata, uk_geom, how="left")
join_left_df

How do I solve this?

I can't seem to make it work.


texts: Open file from QGIS and search for text within it

texts:
I have a lot of data associated with a layer in QGIS. Each individual feature has its own information associated with it - but this information is within a number of different PDF files.

I have managed to create an action that opens the particular PDF file associated with the element easy enough - but I want to be able to search within the PDF once it is open. For example:

open [% "File Link"%]
searchtext = [% "ID" %]

So, I want to open a pdf file from one column (in the attributes table), and then search that pdf with the word in another column of the attributes table.

Essentially, I want the action to "Ctrl F" and jump to the correct page without me having to manually do it.

Alternatively, if I extracted all the PDF's - is there a way to search a folder with the [% "ID" %] name, so I don't have to link every single document path to its corresponding feature?


texts: Why does QGIS load an image icon in layers panel when WFS layer load fails?

texts:
I'm trying to load WFS layers from a GeoServer and having trouble. It works in MapInfo but in QGIS the layer name appears in the layers panel but with an image icon
The layer doesn't load correctly and I've checked the responses from the server using the Network Activity plugin; they all appear to be valid.
What are the reasons why QGIS may load a layer in this way as knowing the answer may help me to understand why the application doesn't seem to work with our particular GeoServer WFS service?


texts: Loading zipped shapefile from URL using PyQGIS

texts:
I am working on importing weather data from an online source into QGIS 3.10 using Python. When I download the data and use the 'vsizip,' it works perfectly fine and loads the file into my project (below).
zip_uri = '/vsizip/C:/Users/camerony/Downloads/day1otlk-shp.zip'
shp = QgsVectorLayer(zip_uri, 'SPC Day 1 Categorical Outlook', 'ogr')
    if not shp.isValid():
        print ("Layer failed to load!")
    else:
        QgsProject.instance().addMapLayer(shp)

However, when I try to use the 'vsizip' and 'vsicurl,' it does not work, and it cannot load the file.
zip_uri = '/vsizip//vsicurl/https://www.spc.noaa.gov/products/outlook/day1otlk-shp.zip'
shp = QgsVectorLayer(zip_uri, 'SPC Day 1 Categorical Outlook', 'ogr')
if not shp.isValid():
    print ("Layer failed to load!")
else:
    QgsProject.instance().addMapLayer(shp)

How would I be able to load this zipped folder and associated shapefiles directly from a website?


texts: Leaflet.js: find the shortest distance between two points (user location and array of points)

texts:
I am trying to work out the shortest distance between the user location and from some collection of coordinates received from an API, and then I used Leaflet Routing Machine to workout the routes. I'm not sure if I am doing it correctly because when I changed the user location, I still get the same route.  Could some one please have a look at my code and advice if I am doing the right thing? I have the code in Jfiddle.

MY APPROACH:

1) When I have received the latlng for user location and from an array of coordinate from an API. I use the Great-circle distance formula Jfiddle
to find the distances between all the points and push them in an array called allDistances

2)After that I look for the shortest distance in the array by using this formula

function closestDistance() {
 return allDistances.reduce((acc, loc) =>
  acc.distance < loc.distance ? acc : loc
 );
}

This is my CODE in  http://jsfiddle.net/jt100/r740b3d2/47/ , please kindly note I have hard coded the user location because all the points are located in the US and my current location is in UK so I can't identify routes between countries.

I am still very new to programming, It would also help me a lot if you could please give me some advice on how I could improve my code. I have been reading about SOLID principle. I want to follow this rule but I’m not sure what I need to change in my code.


texts: Estimating area of overlap between polygons in same layer

texts:
I have a series of buffers representing service areas of bike share stations that extend 400 from the station (network distance). The polygons were generated using Network Analyst. Evidently, buffers of stations that are close to each other overlap. My question is how do I calculate the overlap when the polygons are in the same shapefile (they are rows of an attribute table) and not two different layers?

I tried using the tool 'Polygon Neighbos' but it generated an empty table.

Any thoughts about how I could proceed?


texts:
I am working on Geoprocessing tools and after completing it my tool output is provided as a shapefile with a tool name.

Could you please tell me how to provide the output file with a tool name?


texts: How to feed a raster map with values from other raster map

texts:
I would like to know if there is a tool able to "feed" a current raster map with the values from other raster map (so the values from the current raster map will be replaced by the ones from the other raster map) in order to get exactly the same extent for both raster maps.
I'm using ArcGIS 10.1 for Desktop.


texts: Identifying a co-ordinate system that works in UK?

texts:
I have some data from an old system that contains co-ordinates which I wish to use but I'm unable to establish what sort of co-ordinates they are. They come from the South East region of the UK. For example:

Latitude | Longitude

51181744, 854338

51166014, 1969501

51275331, 707431

51212390, 998914

I have QGIS and ArcGIS available. I've been through a number of tutorials online but have not gotten anywhere with it. It is an American software product but the product is no longer made. I've tried projfinder.com but it isn't working as their server keeps returning a 500 error.

Looking at the extents in ArcGIS shows for all my data:

Top 51467450

Right 2662457

Bottom and left are 0.

I'm aware this is one of the most common questions on here (having already done a lot of reading on this site) but am now out of options. I have a limited understand of this area.


texts: Listing only Polygon Feature Classes in IMxDocument Using ArcObjects?

texts:
I am using VS 2017 with ArcObjects trying to list all polygons in IMxDocument to be used as clip feature in my add-In. I have a WinForms application and the following code:

IMxDocument document = (IMxDocument)ArcMap.Document;
IMap map = document.FocusMap;

for (int index = 0; index < map.LayerCount; index++)
{
    ILayer layer = map.get_Layer(index);
}

I know I will need to use

geometry.GeometryType != esriGeometryType.esriGeometryPolygon

at some point in my code, but I could not manage my code to achieve that.

How can I get only polygons?


texts: Different outputs from the same query when using country vs continent data (OSM and Osmosis)

texts:
I am attempting to extract highways (motorways, trunk, and primary roads) for Europe from OpenStreetMap using Osmosis. As I am new to using Osmosis, as an intial test, I downloaded the country osm.pbf file for Denmark from Geofabrik. Using Osmosis and the follow query, I got the roads I needed:

osmosis --rbf denmark-latest.osm.pbf --tf accept-ways highway=motorway,trunk,primary --used-node --wx roads_denmark.osm

I then downloaded the Europe-latest.osm.pbf from geofabrik and ran the same query expecting the same results, but for Europe. However, there are many roads missing, and when you compare the Denmark output it is clearly incorrect. In the image below, green is the output when the Danish data is used, and red when the Europe data is used. The European data is only a fraction of the output from when the Danish data is used.

My understanding is that if my query works for one dataset, it should work on another. Therefore:

Is there something fundamental I missing in my understanding? Are the country and European geofabrik data different (apart from coverage)?
Are there potential adjustments I can make to my query to extract the highways?


texts: Error while trying to use GRASS GIS' r.reclass module

texts:
I'm not sure I understand how to use r.reclass in GRASS.

Can somebody go through what the inputs should be or direct me to a tutorial (I've looked...) for the current version? I get an error when I try to use it, and I think it's just because my text inputs are wrong or something like that

EDIT:
My inputs were:

Name of raster to be reclassified

grdn36w104_13 [a DEM from the USGS]

File containing reclass rules

(a .txt file with the text)

 1525 thru 3110 = 1
    0 thru 1524 = NULL
    3111 thru 5000  = NULL

Name for raster output map

test

The error message was:

Followed by "Cannot get region map of grdn36w104_13@Elevation_Data", then "Cannot get input region"


texts: google markers are updating but not removing previous markers

texts:
I am working on redo search functionality on map drag event. I am using php as server side language, What I do in this is that when I refresh my code it shows all restaurants in listing and all these restaurants in google map with markers. Now I click on redo search check box and drag map. for this I am using map dragend event, i call ajax function and fetch all restaurants with in the distance of 9 km from middle of the map. I get all records in listing and also in google map as a marker, but the problem is that it does not remove the markers of previous restaurants those are not showing now in listing after search I am also using marker.setMap(null). Here is my code for map.

var map;
var marker;
google.maps.event.addDomListener(window, 'load', initialize);
function initialize(){
    var bounds = new google.maps.LatLngBounds();
    map=new google.maps.Map(document.getElementById("map-canvas"));
    updateMarkers(bounds,map);

    google.maps.event.addDomListener(map,'dragend',function() {
        if($('#redo_search_in_map').is(":checked")){
            formlat = map.getCenter().lat();
            formlng = map.getCenter().lng();
            var type = $("#types").val();
            var keywords=new Array();
            $(".Check").each(function(){
                if($(this).is(":checked")){
                    keywords.push($(this).val());
                }
            });
            $.ajax({
                url:ajax_url+'listing/search_listings',
                data:{
                    'formlat':formlat,
                    'formlng':formlng,
                    'type':type,
                    'keywords':keywords
                },
                type:"get",
                success:function(resp)
                {
                    $("#ajax_listing").html(resp);
                    deleteMarkers();
                    updateMarkers(bounds,map);
                    return false;
                },
                error: function()
                {
                    alert('error');
                }
            });
            return false;
        }
    })

}

function updateMarkers(bounds, map){
    <?php foreach($filterRecords as $key=>$records){ ?>
        var lat = '<?php echo $records['Listing']['latitude']; ?>';
        var longi = '<?php echo $records['Listing']['longitude']; ?>';
        var latlng = new google.maps.LatLng(lat, longi);
        marker = new google.maps.Marker({
            position:new google.maps.LatLng(lat,longi),
        });
        marker.setMap(map);
        bounds.extend(latlng);
    <?php } ?>
    map.fitBounds(bounds);
}

function deleteMarkers(){
    marker.setMap(null);
}


texts: Where in the UK may I find a register or map of properties and their energy use

texts: I am  looking for a register of the energy use of buildings in the uk with granularity down to individual builing. Is this publically available? I recall that this information is collected when processing energy certificates


texts:
I want to convert a coords interval (especially latitude) to real distance.

I mean let's say we stay at the exact same longitude, but with 1 degree of difference in latitude.

How many kilometers/meters/whatever does this 1-dimensional bounding box amount to?

Is there some type to do the conversion?


texts: Processing Problems with the Difference (QGIS) and Erase (ArcGIS) Tool with big data

texts:
I am struggling with processing a huge data set.

I downloaded it from a German geodata site (https://download.geofabrik.de/europe/germany.html). I downloaded the ZIP of Bavaria and clipped the eight layers (landuse, water, traffic, transport, pois, pawf, natural and building) with a vector shapefile which has the  boundary shape of Bavaria. Afterwards I merged all the layers together. The total ZIP has 581MB after clipping and merging the layers the final shapefile has around 4,8GB.

Someone has an idea what happend here?

I used QGIS 2.18.18.

My problem in particular is that I want to get the difference of the clipped and merged layer within the boundary Layer of Bavaria to get a new shapefile with all the empty areas but nothing works. Both GIS programs are calculating for hours and in the end I get an error message. I think the dataset is too big but I also tried to use the difference tool with only the smallest layer and the same happens. So it does work with the huge merged/clipped layer and not with the individual layers.

I worked with the Difference Tool on QGIS and with the Erase Tool on ArcGIS, I also tried the Symmetrical Difference Tool on ArcGIS. The problem is that there is no other possibility for me to get the data of open areas on a different way.

Has someone an idea what could solve my problem?

Is my dataset to big to be proceed with QGIS?

I don't think that my laptop is the problem (I have an Lenovo T470 with 8GB RAM, 64bit QGIS and Core i5 of the 7th Gen) because the computers which I can use at the University with ArcGIS are struggling too.

Here a three pictures as an example.
The first one shows the Bavarian Boundary Layer

The second one shows the clipped an merged Layers

The third picture has the Boundary Layer under the clipped and merged Layer

The yellow dots are my open areas which I need to work with and which I tried to cut out with the Difference, Erase or Symmetrical Difference tools.


texts: Importing coordinates from Google Maps as CSV into QGIS, but points are displaying in wrong location

texts: I'm using QGIS and am trying to import coordinates pulled from Google Maps as CSV but they are displaying in the wrong location. They should be showing up in CA, however they are only displaying in the ocean off the west coast of Africa. I have tried changing projections and CRS but nothing so far has worked.


texts:
I am working on a map of South America (from the top of Chile down to Cape Horn; everything south of 18S).  I need to include distance buffers on my map but I am getting confused by which projection to use to get the distance right.

Would Plate Carree work or is there a better one for this area?


texts: Routing network from database created by osm2pgsql

texts:
I created my PostGIS database by osm2pgsql tool and now I want create network suitable for routing. (I could not use directly osm2pgrouting because I need street names e.g. and so on).

For this network I this is the best solution use pgr_createTopology(). But now I could not use  function for shortest path pgr_dijkstra() because table "planet_osm_lines" not contains column which will be represented cost of edge, which is required for this function.

Can you tell me ho to resolve this?

My table planet_osm_line's script:

CREATE TABLE public.planet_osm_line
(
    osm_id bigint,
    access text COLLATE "default".pg_catalog,
    "addr:housename" text COLLATE "default".pg_catalog,
    "addr:housenumber" text COLLATE "default".pg_catalog,
    "addr:interpolation" text COLLATE "default".pg_catalog,
    admin_level text COLLATE "default".pg_catalog,
    aerialway text COLLATE "default".pg_catalog,
    aeroway text COLLATE "default".pg_catalog,
    amenity text COLLATE "default".pg_catalog,
    area text COLLATE "default".pg_catalog,
    barrier text COLLATE "default".pg_catalog,
    bicycle text COLLATE "default".pg_catalog,
    brand text COLLATE "default".pg_catalog,
    bridge text COLLATE "default".pg_catalog,
    boundary text COLLATE "default".pg_catalog,
    building text COLLATE "default".pg_catalog,
    construction text COLLATE "default".pg_catalog,
    covered text COLLATE "default".pg_catalog,
    culvert text COLLATE "default".pg_catalog,
    cutting text COLLATE "default".pg_catalog,
    denomination text COLLATE "default".pg_catalog,
    disused text COLLATE "default".pg_catalog,
    embankment text COLLATE "default".pg_catalog,
    foot text COLLATE "default".pg_catalog,
    "generator:source" text COLLATE "default".pg_catalog,
    harbour text COLLATE "default".pg_catalog,
    highway text COLLATE "default".pg_catalog,
    historic text COLLATE "default".pg_catalog,
    horse text COLLATE "default".pg_catalog,
    intermittent text COLLATE "default".pg_catalog,
    junction text COLLATE "default".pg_catalog,
    landuse text COLLATE "default".pg_catalog,
    layer text COLLATE "default".pg_catalog,
    leisure text COLLATE "default".pg_catalog,
    lock text COLLATE "default".pg_catalog,
    man_made text COLLATE "default".pg_catalog,
    military text COLLATE "default".pg_catalog,
    motorcar text COLLATE "default".pg_catalog,
    name text COLLATE "default".pg_catalog,
    "natural" text COLLATE "default".pg_catalog,
    office text COLLATE "default".pg_catalog,
    oneway text COLLATE "default".pg_catalog,
    operator text COLLATE "default".pg_catalog,
    place text COLLATE "default".pg_catalog,
    population text COLLATE "default".pg_catalog,
    power text COLLATE "default".pg_catalog,
    power_source text COLLATE "default".pg_catalog,
    public_transport text COLLATE "default".pg_catalog,
    railway text COLLATE "default".pg_catalog,
    ref text COLLATE "default".pg_catalog,
    religion text COLLATE "default".pg_catalog,
    route text COLLATE "default".pg_catalog,
    service text COLLATE "default".pg_catalog,
    shop text COLLATE "default".pg_catalog,
    sport text COLLATE "default".pg_catalog,
    surface text COLLATE "default".pg_catalog,
    toll text COLLATE "default".pg_catalog,
    tourism text COLLATE "default".pg_catalog,
    "tower:type" text COLLATE "default".pg_catalog,
    tracktype text COLLATE "default".pg_catalog,
    tunnel text COLLATE "default".pg_catalog,
    water text COLLATE "default".pg_catalog,
    waterway text COLLATE "default".pg_catalog,
    wetland text COLLATE "default".pg_catalog,
    width text COLLATE "default".pg_catalog,
    wood text COLLATE "default".pg_catalog,
    z_order integer,
    way_area real,
    way "public.geometry"(230633732)
)


texts: Filtering a raster by pixel value ? (Qgis)

texts:
I have downloaded Globcover data in GeoTIFF format. I want to create a new grid which only has values where the Globcover data (the only band in the TIF file) has one certain value (==190).

Can I do this with GDAL?

If not, which way should I go? The procedure should be callable from a Shell script (Linux). I have QGIS installed (but no clue how to use it).


texts: Cartogram not working: "Cancelled by user" - Why?

texts:
I am using the QGIS Cartogram plugin, but unfortunately something is not working as it should.

First, the demo works just fine. But when I want to create a cartogram with my own layer. But as soon as I click „ok“, it says
"Cartogram creation cancelled by user“.

Do you have a clue, what I'm doing wrong? What am I canceling?


texts: Changing name of file download from browser app?

texts:
I am developing browser application using MapServer and JavaScript. I've got part where on map click I get FeatureInfo about an object:

function getFeatureInfoGML(map){
    var openwindowhandle=window.open('', "openwindow")
    if (openwindowhandle){ openwindowhandle.close()}
    map.on('singleclick', function(evt) {
    var viewResolution = /** @type {number} */ (view.getResolution());
    var url = immovableSource.getGetFeatureInfoUrl(
        evt.coordinate, viewResolution, 'EPSG:4326',
        {'INFO_FORMAT': 'gml'});

    if (url) {
                infoWindow=window.open(url, "openwindow", "width=500,height=600");
            }
    });

}

where

    immovableSource = new ol.source.ImageWMS({
        url: 'http://localhost:8000/cgi-bin/mapserv.exe?map=c:/ms4w/project/info.map',
        params: {'LAYERS': 'immovable'},
        ratio: 1,
        serverType: 'mapserver'
        }),

in my map file my output format is:

OUTPUTFORMAT
    NAME "OGRGML"
    DRIVER "OGR/GML"
    MIMETYPE "text/xml; subtype=gml/2.1.2; driver=ogr"
    FORMATOPTION "STORAGE=memory"
    FORMATOPTION "FORM=multipart"
    FORMATOPTION "FILENAME=rezult.gml"
END

When I click on the map a file is downloaded with name "mapserv.exe", but I want it to be *.gml. When I open mapserv.exe it's a GML file, I just want to change the name of the downloaded file.


texts:
SQL help please! I've tried searching, but can't find the help I need.

This is for ArcGIS Advanced 10.2.

I have a parcel layer with all kinds of info that the general public shouldn't be able to access. What I want to do is use a definition query to hide entire columns of data, but still display the parcel with some non-sensitive information. I'd like to avoid deleting the columns if I can.


texts: Filter geoserver layers based on postgis geometry

texts:
I am wondering how, and if it is possible, to filter a geoserver's layer (from a large shapefile) on a specific geometry coming from a postgis table.

I have a huge shapefile for a country. I'd like the user to only be able to view a limited area of the layer he has been granted permission for (a city or whatever). The geometry he is allowed to see is stored as a multipolygon in a postgis table.

Is there a way to do that directly in geoserver (server-side)?

If not, could you give me some hints / technologies to use?