texts: Replace few symbols in all strings using QGIS Field Calculator?

texts:
I created the "Photos" column in the attribute table. Each line contains a unique path to a specific photo.
For example:

D:/.../IMG0031.JPG

D:/.../IMG0048.JPG

Now I need to replace the path for each photo from drive D to drive C.

C:/.../IMG0031.JPG

C:/.../IMG0048.JPG

At the same time, the right side in each line should not change!
What function should I use?


texts: How to export vectors from QGIS to MapInfo?

texts:
I am trying to save some vector lines that I made in QGIS (1.8 Lisboa) to a format that can be added to MapInfo Professional 6.0 (ancient version I know..).

When I save the vector to a Mapinfo format it produces two files: .MIF & .MIDI.

However in the Mapinfo if appears as though a vector file has 5 separate files:
.ID .IND .MAP .MIF .DBF......and sometimes a .MIDI file.

Can I not add the vector line without these other files? If so is there any way that I can produce them? Thanks!


texts: How to generate smallest possible country boundaries file

texts:
The problem i'm trying to solve is rather simple:
Given any point X lookup the country Y which it belongs to
Assumptions:

Country boundaries are expected to be intersection free and have a reasonably high resolution
All points being looked up are part of a landmass and should therefore be assignable to a country

I'm currently using the admin0 dataset provided by the NaturalEarth project to solve this but wondered if it can be simplified further to save space and computation time without losing too much precision.
My first naive approach was to apply a line simplification algorithm using https://mapshaper.org/. While this was successful to a certain degree, i noticed a negative impact on the resolution and landmass coverage:

Douglas-Peucker

Visvalingam

Would it be possible to only simplify non-shared line segments like coastlines and group islands? Geofabrik seems to be doing something similar when displaying the extents of their downloads:


texts: Any good ruby gems that can calculate a polygonal geofence?

texts: Looking for a Gem that has polygonal geofencing. Many gems offer square or radius geofencing.


texts:
I have a script that runs Tabulate Area on a single feature class multiple times with different rasters.  In total Tabulate Area runs 5 times.  2 of the input rasters are provided as arguments to the script and Tabulate Area works as expected.  The other 3 rasters are results of raster functions (i.e. Reclassify, SetNull, EucDistance, etc.).

Every time I run the script, Python crashes when trying to run Tabulate Area on one of those preprocessed rasters.  It happens on a different raster every time, seemingly randomly.  There is no error message, Python simply crashes.  When I hit Get Details in the generic Windows error dialog, it seems like there is a different Fault Module every time as well.

When I try to run Tabulate Area manually in ArcCatalog, ArcCatalog crashes and asks if I want to send a report to ESRI.  Again, no error message.

UPDATE:

Here is some of the code involved:

def GetBuffer(inRaster, lcClass, bufDist, ringBuffer=True):

    setnull = SetNull(inRaster, 1, "Value <> " + str(lcClass))

    dist = EucDistance(setnull, bufDist)

    buf = Reclassify(dist, "Value", RemapRange([[0, 100, 1], ["NODATA", 0]]))

    if ringBuffer:

        reclass = Reclassify(setnull, "Value", RemapValue([[0, 0], [1, 1], ["NODATA", 0]]))

        ring = arcpy.gp.Minus_sa(buf, reclass)

        return ring

    else: return buf

def SubtractRaster(inRaster1, inRaster2):

    sub = arcpy.gp.Minus_sa(inRaster1, inRaster2)

    setnull = SetNull(sub, 1, "Value <> 1")

    return setnull

lcRaster = arcpy.Raster(landcover)

treecanopy = GetBuffer(lcRaster, TreeCanopyClass, TreeCanopyBuffer, False)

parkinglots = SubtractRaster(GetBuffer(lcRaster, ParkingLotClass, ParkingLotBuffer), treecanopy)

buildings = SubtractRaster(GetBuffer(lcRaster, BuildingClass, BuildingBuffer), treecanopy)

impervious = SubtractRaster(GetBuffer(lcRaster, ImperviousClass, ImperviousBuffer), treecanopy)

...

TabulateArea(fc, "GID", parkinglots, "Value", plTabArea , 1)

TabulateArea(fc, "GID", buildings, "Value", buildTabArea , 1)

TabulateArea(fc, "GID", impervious, "Value", impTabArea , 1)

UPDATE 2:

I installed ArcGIS 10.1 Service Pack 1 after seeing bugs related to Tabulate Area in 10.1.  When I first ran the script with SP1, it worked.  When I ran it again later, I had the same issue where Python simply crashed without providing an error message.


texts: OSM update without losing own additional data

texts: I have added data (e.g. name keys) to my OSM PostgreSql database. How can I update my database without deleting or overwriting this data? In other words, how can I merge newer OSM data to my existing database without losing my additional data?


texts: Issue rendering polygons at global scale when using the Robinson coordinate reference system in QGIS 2.10

texts:
I am trying to create a map using the Robinson coordinate reference system (ESRI:54030) in QGIS 2.10, on Mac. I used the 'on the fly projection' feature in the Project Properties.

My base layer is a polygon vector file of the world's countries.
Some of my other vector layers render correctly at any scale (points, lines), but the world countries shapefile does not render correctly when visualizing the map at its full extent. A bounding box is rendered instead of the country's boundary.

Here is a picture of the issue (full extent):

As I zoom in, the polygons render correctly. I've tried with various world boundaries datasets from various sources but the issue remains.

Anyone has encountered this problem before ?


texts: Loading a raster using Openlayers and Mapserver

texts:
I'm trying to load a raster into my browser using Openlayers. I can examine my WMS raster layer using OpenJump, but I'm not able to get it to load into my browser. My mapfile is as follows:

MAP
    ###########
    NAME         "ny_pop"
    SHAPEPATH     "/v0.3/mapserver/"
    UNITS         meters
    EXTENT        -8269337 4955373 -8190365 5005745
    SIZE          1000 1000
    STATUS        ON
    DEBUG         ON
    IMAGETYPE     "png"

    ###########
    PROJECTION
        "init=epsg:900913"
    END

    ###########
OUTPUTFORMAT
  NAME "png"
  DRIVER AGG/PNG
  MIMETYPE "image/png"
  IMAGEMODE RGB
  EXTENSION "png"
  FORMATOPTION "GAMMA=0.75"
END

    ###########
    WEB
        METADATA
            "wms_title"           "WMS Pop Tif"
            "wms_onlineresource"  "http://mysite.com/mapserv.cgi?map=raster_1.map&"
            "wms_enable_request"  "GetMap GetFeatureInfo GetCapabilities"
        END
    END

    ###########
    LAYER
        NAME      "ny_pop"
        DATA      /home/mysite/webapps/htdocs/v0.3/mapserver/ny_pop.tif
        STATUS    ON
        TYPE      RASTER
        METADATA
            "wms_title"  "WMS NY Raster"
            "wms_srs"    "EPSG:900913"
        END
    END

    ###########
END

Everything (map, raster and mapfile) are using EPSG:900913, and I have an existing map (centered around New York) which I am adding the layer rasterLayer to. My javascript is as follows:

function generateRaster(){

    console.log('get map bounds');
    bounds = map.getExtent();

    rasterLayer = new OpenLayers.Layer.WMS(
           "ny_pop",
           "http://mysite.com/mapserv.cgi?",
           {map: 'raster_1.map'},
           {'maxExtent': bounds}
       );
    console.log(rasterLayer);
    map.addLayer(rasterLayer);

}

I'm trying to limit the raster to the bounds of the map using the {'maxExtent': bounds} command, but I also tried with just the map name and URL and it did not work either. I'm still trying to get my head around the .map format, so there could be some mistakes there also. However, as I could view the raster in OpenJump I assumed it was ok.

EDIT:
I updated my javascript to the following:

rasterLayer = new OpenLayers.Layer.WMS(
   "NY Population layer..... ",
   "http://urbmet.webfactional.com/mapserv.cgi?",
   {map: 'raster_1.map', layers: 'ny_pop'}
);

but still nothing shows up on the map. I've added a screencapture of my firebug output - this is the contents of console.log(rasterLayer).


texts: Copying Personal Geodatabase (in Access) with relationship classes to ArcSDE for SQL Server?

texts:
How to copy a Personal Geodatabase (in Access) with relationship classes to ArcSDE for SQL Server.

I tried a simple copy paste and imported, but didn't get the relationship classes.

  error: unable to initialize weights in feeder manager.

ArcGIS Server 10.0 (ArcSDE), SQL Server 2008 r2


texts: Why is st_write not saving this dataframe as a shapefile?

texts:
I have been working with sf package for several weeks now. I have been able to figure out the solution of problems by myself or searching on forums. However, this one has been demanding too much time and I don't know what is the problem yet.
The task is quite basic, I want to write a dataframe as .shp with the following code, which is pretty straightforward:
st_write(UPS, "Shapefiles/UPS.shp", layer_options = "GEOMETRY=AS_XY", delete_layer = TRUE)

Some of the variables are the following (only SHAPE_Leng, SHAPE_Area and geometry are geographical).
And these are the descriptions of the variables (I cannot generate a random database to replicate since this specific dataframe has been quite a trouble):

When I execute the aforementioned code I get the following output:
Deleting layer `UPZ2' using driver `ESRI Shapefile'
Writing layer `UPZ2' to data source `Shapefiles/UPZ2.shp' using driver `ESRI Shapefile'
options:        GEOMETRY=AS_XY
Updating existing layer UPZ2
Writing 111 features with 11 fields and geometry type Multi Polygon.
Unknown field name `UPlCd': updating a layer with improper field name(s)?
Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :
  Write error
In addition: Warning message:
In abbreviate_shapefile_names(obj) :
  Field names abbreviated for ESRI Shapefile driver

So, as you can see, I don't have any variable named UPlCd, so I don't have figured out yet how to proceed. In addition, I don't know what is Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :  is referring to since no variable is named like that.
The question is: how can I succesfully write this dataframe as a .shp file? (and in the process, which details I'm missing).
I have seen these responses (1, 2 and 3) but they have not been helpful to clarify which problem I'm having with this dataframe so I cannot save it.


texts: Angle offset in the Line Direction Histogram Plugin

texts: Can someone please explain what the angle offset is in the QGIS Line Direction Histogram Plugin. I have tried using a search engine to figure it out but I cant. I am using the plugin to calculate orientation of different vector lines.


texts:
I am using OL3 to load a geojson layer with points on a map. Although the layer is loaded and it can be seen, when I zoom out the position of the points change on the map. I thought this would related to the projections of the map and the points but I can not figure out how to fix it.
This is the script:

var format = new ol.format.GeoJSON();
var url = 'https://gist.githubusercontent.com/cmiles74/ef62e43595deabb68f5a/raw/ef85f8e7ed34ccbf5768ade25e4d3485ded16938/gistfile1.json';

var vectorSource = new ol.source.Vector({
strategy: ol.loadingstrategy.bbox,
loader: function(extent, resolution, projection) {
  $.ajax(url).then(function(response) {
      var features = format.readFeatures(response,
        {featureProjection: 'EPSG:3857'});
      vectorSource.addFeatures(features);
  });
}
});

var map = new ol.Map({
layers: [
  new ol.layer.Tile({
      source: new ol.source.XYZ({
        url: 'http://server.arcgisonline.com/ArcGIS/rest/services/' +
            'World_Topo_Map/MapServer/tile/{z}/{y}/{x}'
      })
    }),
  new ol.layer.Vector({
    source: vectorSource
  })
],
renderer: 'canvas',
target: 'map',
view: new ol.View({
  center: [-9343811.351360613, 5424617.381580625],
  zoom: 14
})
});

The geojson looks like this:

And in the below images you can see how the position of the points change in different zoom levels:


texts: Prepared statements and MULTILINE geographies

texts:
I'm working on a web application that parses GPX files and enters the contents into a database as Geography objects. Seeing how incoming files can be potentially insecure, I was hoping to use prepared statements for this. These files can contain multiple track segments, so I am using MULTILINESTRING as the content type.

Normally I would make a geography object in db (postgres) using something like this:

ST_GeographyFromText('MULTILINESTRING((lon lat, lon lat), (lon lat, lon lat, lon lat))')

I will create the sequence of lines and points by iterating through the GPX object, but what's a good way to make the insertion statement safe? I tried:

$result = pg_prepare("query", 'insert into geogs (title, setID, geog, order) values ($1, $2, ST_GeographyFromText($3), $4) returning id');

But the $result is a boolean false. Is there any way to use prepared statements to achieve this?

EDIT: that was one silly code bug. The problem was not at all related to GIS or prepared statements. The db was just upset that I used order as a column name. Otherwise the whole this works just fine after renaming it:

$result = pg_prepare("query", 'insert into geogs (title, setID, geog, sequence) values ($1, $2, ST_GeographyFromText($3), $4) returning id');


texts: How to improve rendering performance of a very large shapefile?

texts:
I have a polygon shapefile that has a 100MB .dbf, and 500MB .shp file component. The reason it is so large, is that it is  classified as a landbase for a whole district.

Everytime I view the file in ArcCatalog, or ArcMap, and move the view window slightly, the entire file needs to redraw from scratch. I have tried spatial indexing and importing into a geodatabase - neither approach provides any noticable performance improvement with respect to rendering.

Esri's help page suggests that that to improve shapefile performance, the user can generalize the file. While this would obviously work, I do not want to lose any information. Splitting the file up is not ideal since I am doing a lot of geoprocessing/queries with its entire area. I guess I could avoid viewing the entire area at once - but sometimes, for example, it is good to see what parts of the file a query has selected.

Is there any other approach I could take to improve rendering performance?

(In theory, building shapefile "pyramids" would be ideal - I'm unsure why ArcGIS has never supported such an approach  - at least that I am aware of...)


texts: Distance in meters between two geographic coordinates

texts:
I want to determine the distance in meters between two geographic coordinates. The points are:

     Lat (deg)   Lon (deg)
A: -33.3930906 -70.5452765
B: -33.3993212 -70.5468767

So if i calc the distance in google maps or with a python library the measure is approx. 700 m but when I do the same query (distance) in PostGIS it tells me it is 291 (I don't know the unit, I believe that is meters)

My queries are:

SELECT ST_Distance_sphere(
    st_makepoint(-33.3930906, -70.5452765),
    st_makepoint(-33.3993212, -70.5468767));

SELECT ST_Distance_Sphere(
    ST_Centroid(ST_GeomFromText('POINT(-33.3930906 -70.5452765)',4326)),
    ST_GeomFromText('POINT(-33.3993212 -70.5468767)',4326));

SELECT ST_Distance(gg1, gg2) As spheroid_dist,
       ST_Distance(gg1, gg2, false) As sphere_dist
    FROM (SELECT
        ST_GeogFromText('SRID=4326;POINT(-33.3930906 -70.5452765)') As gg1,
        ST_GeogFromText('SRID=4326;POINT(-33.3993212 -70.5468767)') As gg2
    ) As foo  ;

And each of the three queries give me a result of approx. 291 or 292.

So, how can I do it?

My original intention is to query all the points in a certain radius. For this I'm using ST_DWithin.

And in my database I store it like:

// Column is
coordinates GEOGRAPHY(Point, 4326)  NOT NULL,
// Insert is
ST_GeographyFromText('POINT(-33.3940301 -70.5482872)', 4326)


texts: Connectors do not appear automatically in easy custom labelling

texts: I'm using the easy custom labelling plugin in QGIS 3. When I move a label away from its associated polygon a connector does not automatically appear as in QGIS 2.18. How do I get connectors to appear in QGIS 3?


texts: Using PyQGIS to append index to attribute table

texts:
How can I add a new column to an attribute table that just numbers each row 1-n with Python in QGIS?
Normally I just grab the code from the log history in QGIS, but there's no output for field calculator.
Better yet, does the attribute table an invisible, implicit index, such that I can access row numbers without having to add a field for it in the first place?


texts: Running GDAL from a docker container using Python

texts:
Dockerfile

FROM python:3.6-stretch

RUN pip install --upgrade pip

# Add unstable repo to allow us to access latest GDAL builds
RUN echo deb http://ftp.uk.debian.org/debian unstable main contrib non-free >> /etc/apt/sources.list
RUN apt-get update

# Existing binutils causes a dependency conflict, correct version will be installed when GDAL gets intalled
RUN apt-get remove -y binutils

# Install GDAL dependencies
RUN apt-get -t unstable install -y libgdal-dev g++

# Update C env vars so compiler can find gdal
ENV CPLUS_INCLUDE_PATH=/usr/include/gdal
ENV C_INCLUDE_PATH=/usr/include/gdal

RUN pip install \
    GDAL==3.0.4

RUN mkdir -p /app/

COPY / /app/

WORKDIR /app/

Running a python script with:

def colour_relief():
    subprocess.call(['gdaldem', 'color-relief', dem_imput,
        r'/app/color_relief.txt',
        colour_relief_output,'-b','1' ])

Results in the below error:
FileNotFoundError: [Errno 2] No such file or directory: 'gdaldem': 'gdaldem'

I'm unable to get any results from gdal --version from the bash command line.

I am looking for the simplest way to run python with gdal from a container.

Original dockerfile is from:
https://github.com/thinkWhere/GDAL-Docker/blob/develop/3.6/Dockerfile


texts: Get lat lon extent of country from name using Python?

texts:
Is there any open source GIS function in python to get lat lon window given only country name. Something like:

get_country_window('United states of america')

should return:

[-135, -50, 22, 50]


texts: Using ArcGIS and QGIS together in multi-user editing environment?

texts:
We have a GIS staff of 7, and another group of about 14 who "use GIS" from time to time. We really, really want multi-user editing. Mostly of a large point file, and we're almost never in the same location, so editing the same point is not a concern.

I'd prefer everyone on ArcGIS, using ArcSDE, but we don't have the budget, or using QGIS and ArcGIS on PostgreSQL/PostGIS, but ArcGIS for Desktop cannot edit PostGIS without ArcSDE.

Currently, I've got the GIS staff on ArcGIS Advanced, and the others on QGIS 2.2. It's working well, multiple QGIS users are able to edit at the same time + one additional Arc user, but I'm always concerned about something popping up that I haven't considered. I'm backing up our data 4x/day.

Hoping there may be others in a similar environment who could help me learn from their experience, avoid making mistakes or losing data.


texts: Permission Denied when installing openstreetmap server

texts:
I am installing a tile server from this link on my Ubuntu Machine 14.04. The problem I am facing is in this line :

./make.py
cd ../OSMBright/
carto project.mml > OSMBright.xml

Whenever I enter 'carto project.mml > OSMBright.xml' it is giving me permission denied error. I have also tried with sudo carto project.mml > OSMBright.xml but still it is giving me permission denied. I am stuck here at this stage and I need to install this server. See I have tried all combination to execute this command. Could anyone help with this?

I tried to post this question on Ubuntu forum but they don't have the proper tags that I am looking for.


texts:
QGIS has been very good about prompting me to update. Not this time. I run QGIS 2.18.12 on Win10. From reading this board I realize that 2.18.13 has been out for at least a week. Not only does QGIS (which I run daily) not prompt me to update, but the Version Information panel states that 2.18.12 is the current version ...

... whereas qgis.org says otherwise. Why is that?

[UPDATE] Coincidentally (or not), about an hour after my original post, the issue resolved itself -- i.e., the Version Information panel knew that 2.18.13 was available, and the yellow update notification message appeared in the status bar. Still unsure what happened and why, so I'll leave this open for now. Maybe someone will offer insight into this mystery.


texts: Show Log of the Sql statements that QGIS generates to retrieve Data from Postgis

texts:
Is there a ways to
- make QGis log its generated SQL Statement
or to
- force POSTGIS to log all interaction with QGIS?

There is something strange happening to the WFS layer I retrieve in QGIS, since it does  not select all data available in the respective database table. I need to find out, what criteria QGIS chooses in order to narrow down the select set.


texts:
I have a plugin which finds buffer for list of cities, provided the distance. The state and city names are taken from the attribute table and gets filtered accordingly. What I want is My plugin should identify the layer name or order of the layer in canvas,irrespective of other layers present in the canvas and access the corresponding attributes from that layer.

I am also just curious whether pointing a specific layer name in code would cause any error in iteration though some other layers are present?

Below is my code please tell me where should I make changes and what would be the change?

    if dist and centerCity:
        #QMessageBox.information(self.dlg, "info", "both True")
        st = '"name" = \'' + centerCity + '\''
        exp = QgsExpression(st)
    else:
        QMessageBox.warning(self.dlg, "Enter the distance","Enter the distance and try again.")
        return          #terminate the function

    layer = self.iface.activeLayer()
    it = layer.getFeatures(QgsFeatureRequest(exp))
    feature = it.next()
    mbuf = feature.geometry().buffer(dist, 2)

    iterFeat = layer.getFeatures()

    for f in iterFeat:
        geom2 = f.geometry()
        valTest = QgsGeometry.within(geom2, mbuf)


texts: I downloaded a bunch of wpts from my gps. I converted them to .SHP format with a CRS WGS 84  in QGIS, and started a new project. I then exported some satelite imagery from the internet as a tiff and imported it as a raster in my project as WGS 84. I have CRS 'on the fly, enabled, but I still cannot see both layers at once. I can zoom between layers but they wont overlay.


texts: While developing stage, Is Proxy.cgi required to get Openlayers WMSGetFeatureInfo?

texts:
I have published the data in Geoserver (Tomcat and OS Ubuntu) and I am developing a map application in my local machine (OS windows xp) using asp.net. In my application, I can get WMS feature information in IE8 with the confirmation msg "This page is accessing information that is not under its control. This poses a security risk. Do you want to continue?" given Yes. But in Chrome, i cannt get feature info. It says

XMLHttpRequest cannot load http://xxx.xxx.x.xxx:8080/geoserver/cite/wms?LAYERS=cite%3Aplym_p1_examplec…e%2Fpng&INFO_FORMAT=application%2Fvnd.ogc.gml&SRS=EPSG%3A27700&X=596&Y=160. Origin http://localhost:3890 is not allowed by Access-Control-Allow-Origin.

I have found similar questions Using WFS protocol in OpenLayers? but its for after hosting the web application in geoserver.My case is still its in under developing in my local machine. My code is:

roadidPicker = new OpenLayers.Control.WMSGetFeatureInfo({
url: 'http://xxx.xxx.x.xxx:8080/geoserver/cite/wms',
    title: 'identify features on click',
    queryVisible: true
});
roadidPicker.infoFormat = 'application/vnd.ogc.gml';
roadidPicker.events.register("getfeatureinfo", this, pickRoadid);
osMap.addControl(roadidPicker);
roadidPicker.activate();

function pickRoadid(e) {
if (e.features.length > 0) {
    var val = e.features[0].attributes['theme'];
    alert(val);
}
}

I have configured the proxy.cgi file in Geoserver (usr/lip/cgi-bin/proxy.cgi) and added my local ip address in

allowedHosts = ['xxx.xxx.x.xxx', 'www.openlayers.org', 'openlayers.org', 'labs.metacarta.com', 'world.freemap.in', 'prototype.openmnnd.org', 'geo.openplans.org', 'sigma.openplans.org', 'demo.opengeo.org', 'www.openstreetmap.org', 'sample.azavea.com', 'v2.suite.opengeo.org', 'v-swe.uni-muenster.de:8080', 'vmap0.tiles.osgeo.org','www.openrouteservice.org']

and added the below line in my script.

OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

but it says another error in chrome console:

GET localhost:3890/usr/lib/cgi-bin/….…T%3Dapplication%252Fvnd.ogc.gml%26SRS%3DEPSG%253A27700%26X%3D682%26Y%3D284 404 (Not Found)

I dont know whether its right or wrong. Proxy.cgi should be copied in Geoserver or my machine? in which path?

Please anyone suggest me to fix this issue


texts: Network Dataset connecting link file - no links

texts: I've produced a network from ITN and Urban Paths data, and both sets of data work separately, however you cannot route between the paths and the road network because no connecting links were created. Has anyone had this problem before?


texts: Creating Thick Neatline using ArcGIS Desktop?

texts:
How would I change the thickness of a neatline using ArcGIS 10.2.1?

Alternatively, if neatline thickness can't be changed, how would I draw it with a rectangle with no fill and rounding?

The map view is a rectangle occupying all but the top and the bottom of the neatline which have the title and legend respectively.


texts:
I have two GeoDataframes - one has polygons, the other points.

import geopandas as gpd
from shapely.geometry import Point
area = gpd.GeoDataFrame.from_file('areas.shp')
print area.crs

{u'lon_0': -91.8666666667, u'datum': u'NAD83', u'y_0': 3000000, u'no_defs': True, u'proj': u'lcc', u'x_0': 6200000, u'units': u'm', u'lat_2': 77, u'lat_1': 49, u'lat_0': 63.390675}

Addresses were assigned geolocations through an API, so I assigned a CRS 4326, as that's the one most commonly used it seems.

location['geometry'] = location.Geolocation.apply(lambda x: Point(eval(x)))
gdf = gpd.GeoDataFrame(location, geometry='geometry', crs={'init': 'epsg:4326'})

Then I changed the area's CRS:

area = area.to_crs({'init': 'epsg:4326'})

And tried to merge the two dataframes

merged = gpd.sjoin(area, gdf, op='intersects')

But the result it gives me is Empty GeoDataFrame. I am not sure why, because the points definitely do fall inside polygons. What am I doing wrong?

Here are my files

https://www.dropbox.com/sh/9lioawgsor5y4wc/AABkaaFpjAjisQEVaLrFwSnsa?dl=0

I suspect that its because of the points that were somehow not properly processed into the right coordinate system maybe.


texts:
If I wanted to create an add-in for ArcGIS that exported the styling of a layer to SLD how much work would it be and what languages could I use?

More info:

I'm an experienced programmer, and have done a lot of coding for GIS systems but am not a user of ESRI products and have never tried writing an extension for any of their products.

I would only need to export basic styling to SLD, not write a comprehensive tool that works for everything.


texts: Generating random valued random points using ArcPy?

texts: I want my area filled radomly with a number of  points with a number of random numerical values in the range i specify. Other than tediously creating each point and filling up random values, is there any better method?


texts: Merge polygons based on five first digits in id

texts:
I have a shapefile with 1364 polygons in the Stockholm area:

But the "areas" are too detailed - I want to make fewer and bigger polygons. I would like to merge/dissolve all "BASKOD2010" id's based on the five first digits. As an example from my table above, I want to merge 3510280, 3510270 and 3510260 into one polygon.

Can anyone explain how I do this in QGIS?


texts: Create binary raster with a spatial condition for generating cell values

texts:
I have a raster grid. Cells either have the value 0 (blank in the example image below*) or NA

I need to designate some non-NA cells with the value of 1 with the following rule:

Cells with a value of 1 cannot touch another cell with a value of 1 (i.e. the surrounding 8 cells must have a value of 0 or NA)

Here is a visual example, where cells with a value of 1 are black

My problem is I have no idea how to generate this using ArcGIS or R. Does anyone have any ideas on how to do this?
*The actual grid has ~15000 cells


texts: How to randomize colors by class in QGIS new symbology?

texts: As a new user of QGIS, one of the tasks I frequently perform is creating a layer where polygons are colored by class, classified by some nominal field in the attribute table. I have developed a habit of switching to the old symbology, because it makes this task easy for an arbitrary number of classes and I can re-randomize the colors until I'm happy with the results. I haven't found any way to do the same thing with the new symbology. The closest I've found is to select "categorized" in the new symbology and select a color ramp, which works well for ordinal data but not so much for nominal categories. I presume that such an important function was not lost in the design of the new symbology, but that I've simply been unable to find it in the manual. I'd prefer to work with the new symbology, simply because that's the direction that QGIS development is going. Can someone explain or direct me to appropriate documentation of how to accomplish this? Thanks in advance.


texts:
I've looked around the net and I've looked for examples, but I just can't seem to find a way to make a WFS layer be the baselayer. I set up the projection the right way and I set isBaseLayer: true . When I open my map it's empty. Only the controls are present and there is a white map.

Could it be that it's simply not supported (but since there is the option to set it as isBaseLayer: true, this wouldn't make no sense) ?!

Thanks for your help.

P.S.: The only thing I've changed can be seen inside of the code.

EDIT:

Here is my code:

function myInit2(){

    var geographic = new OpenLayers.Projection("EPSG:4326");
    var mercator = new OpenLayers.Projection("EPSG:900913");
    var world = new OpenLayers.Bounds(-180, -89, 180, 89).transform(
            geographic, mercator
    );
    var deutschl = new OpenLayers.Bounds(5, 50, 15, 56).transform(
            geographic, mercator
    );

    var dortmund= new OpenLayers.LonLat(7.46, 51.51).transform(
            geographic, mercator
    );

    map = new OpenLayers.Map("map_location", {
        maxExtent: new OpenLayers.Bounds(5, 50, 15, 56).transform(
                geographic, mercator)
    });

    //var osm = new OpenLayers.Layer.OSM();
    // this was the old baselayer

    var baseL = new OpenLayers.Layer.WFS(
                "Bundesländer",
                "http://vs136213:8080/geoserver/wfs",
                {typename: "BitmarckKarten:land_shp"},
                {isBaseLayer: true, extractAttributes: true, styleMap: myStyleMapL}
        );

    //initialize the base layer here
    //map.addLayers([osm]);
    //old add of baselayer

    map.addLayers([baseL]);

    map.addControl(new OpenLayers.Control.LayerSwitcher());
    // map.addControl(new OpenLayers.Control.Permalink());
    map.addControl(new OpenLayers.Control.MousePosition());
    map.addControl(new OpenLayers.Control.PanZoomBar());
    map.addControl(new OpenLayers.Control.Navigation());
    map.setCenter(dortmund, 6);
    getDivInfos();

};


texts:
I have the following query that finds all id inside POLYGON(({points}))

SELECT id
FROM t
WHERE ST_DWithin('POLYGON(({points}))', ST_Point(latitude, longitude), 0);

The query is slow because it doesn't use the (latitude, longitude) index and must compute the formula for every possible pair of points.

How can I change my query to force Postgres to use the (latitude, longitude) index (I need to change the query because I cannot add some other index or create/change the table)?

I have the following index:

"index_latitude_longitude" btree (latitude, longitude)

latitude, longitude has double precision type.

I think that Postgres will use the index if we add something like: latitude <= ... and longitude <= .... How can we do this?


texts:
My leaflet map has two WMS layers with different scales; one of them is a 4cm topographic map. I would like to have both WMS on the map but showing only the topographic map when the zoom of the map is sufficient to read something.
Currently I use the minZoom / maxZoom attributes in order to switch when zoom 8 is reached:

var topo4cm =L.tileLayer.wms(topoaddress, {
    service: "WMS", transparent: !1, minZoom: 8, maxZoom: 18
});
var carto=L.tileLayer(cartoaddress, {
        minZoom: 0, maxZoom: 7, zoom:function() {
            var e=map.getZoom();
            return 10>e?"L0"+e: "L"+e
        }
    });

This does what I want: when zoom = 8, the topo map is shown.
However the other layer (carto) cannot be manually switched in the panel control.
Would there be a way to allow the user to manually switch back to carto?


texts: Extracting very small local ridges out of DEM

texts:
I have a DEM with small pits (one pit is approximately 5-10 cm wide, here shown as Hillshade).

I want to extract the ridges of these pits to have polygons out of each pit (like this where I did it manually)

The main problem is, that the for one pit the ridge on the right side f.e. is much higher than the one at the left side, because the overall of the wall is quite curved and seen from the front then the ridge on the left could have the same value than the "deepest" point of the pit. I tried to visualize what I mean - that is the wall drawn if you see it from the top.

So, on the one hand I tried to calculate the bigger mean of the whole area and then to nivellate the whole DEM to get it like in the picture above on the right side. However the problem is that my whole area is not only going convex but also concave and in all possible directions so that it is not possible to nivellate it properly.

The solution I reached up to now is, that I calculate the curvature and then thinning the areas. However, the thinning resolved on lines being not properly on the ridges and some weird polygons that do not match at all - see below. I also tried it with Topographic Position Index, however, the ridges are not marked well enough in many cases to close the pits into single polygons.

Does someone have another idea how to be able to extract those ridges as good as possible?

I am open to all kinds of software and solutions: ArcGIS, QGIS, R, Python,...


texts: QGIS 3x automeasure distance for all lines in the layer

texts:
I have got 2 lines in my layer. I would like to measure their distance automatically.
I found the tool, which might be helpful. This is the v.distance algorithm in GRASS, as advised in the link below:
QGIS calculate distance of point along a line
My data attribute table looks like this by far:

and I would like to have the distance for all of them in the Length [meters] column.
I tried the v.distance GRASS plugin, but the situation looks as follows:

I am getting an error:
Wrong or missing parameter value: Column name(s) where values specified by
'upload' option will be uploaded

and I don't really know what to do with it.
I found some similar problem, where the v.distance algorithm also didn't work properly,
GRASS v.distance in QGIS not working
but they suggest changing the code in Python section. Is it still needed?
Is there any way to autocalculate the distance for all line features available in the given layer?


texts: Convert polygons in shapefile to a NumPy array?

texts:
I have a shapefile containing a number of polygons and all attributes being text. I'm trying to import the polygons into NumPy as an array where each polygon is represented as unique values.

I approach this by using gdal_rasterize to generate a GeoTIFF, which I then can convert to an array:

gdal_rasterize -a provcode -l longhurst PG:'host=localhost dbname=biomes' -tr 0.25 0.25 out.tif

and

tif = gdal.Open('out.tif')
tifArray = tif.ReadAsArray()

My problem is that all polygons get the same values since provcode is of type string. How can I make gdal_rasterize burn different values to represent different polygons?

Alternatively, is there a better way to do this conversion?


texts: How to do a Kernel Density map from weighted point observations using a Cauchy kernel?

texts:
I need to produce a kernel density map, but instead of using a common kernel function such as Gaussian, I need to use the Cauchy kernel function.

I have tried the Kernel Density tool in ArcGIS and Heatmap plugin in QGIS, but neither of them have the option to use the Cauchy kernel function.

I am new in GIS and R, and I don’t know how to proceed.

Please let me know how this can be done.


texts: Label only specific features related to Atlas feature?

texts:
I have a shapefile of about 100 streets and another one of all the districts within a city. Some of these streets run over multiple districts. In QGIS, I use the atlas tool to successively print maps of theses streets one by one.

How can I only display the labels for the districts under each individual street street in the atlas?

I tried doing this with rule based styling. I tried using the refFunctions pluging and other spatial operators in the field calculator.


texts: ArcGIS Server Services status 'stopping'

texts: My services show as 'stopping' in the arcgis server 10.3.1 manager after the upgrade from 10.2. In the admin page the status reads as started for the configuration state and stopped for the real time state. I can't publish a new service as there is error that my publishing service isn't started and i can't start it using any of the available ways. Is there anyone with the same issue?


texts: Enabling an Edit on Shapefile in QGIS only with a password?

texts:
I´m customizing some tasks in QGIS 2.14 using Python.

I am managing these edits in different layers (always shapefile layers, not database).

In one of them I want the users to enter a password before they can edit.

How can I implement this?


texts: Conflicts using Tkinter/wxPython with Python AddIn scripting?

texts: I want to design a custom dialog window where i can select various datasets of choice using scroll down option and specify some conditions and thresholds. Right now i am using python addins within arcgis 10.1 for toolbar and tools. I came across tkinter and wxpython as options for dialog box designing. Will there be a conflict between using tkinter or wxpython within addin scripting?


texts: Calculating distance point to point with certain criteria based on field using ArcMap?

texts:
I have a dataset of points representing solar parks, also these points contain data on the year the solar park was developed.
Now I would like to know the distance of these solar parks to the nearest solar park at the time it was developed.

I know how to get the distance between nearest points, but I would like to add a criteria based on the field "year developed". So that I dont get the distance to the nearest point, but to the nearest point with a value for "year developed" lower (earlier in time) than the point the distance is calculated from.

In short, I need to calculate distance to nearest point with a value for "year developed" < than value for designated point.

Does anyone know how to do this using ArcMap 10.3?


texts: I have a lot of numbered classes in my eCognition project (version 9.0.1), which contains approx. 2 Mio. segments. All segments have been assigned to a class and now I want to merge them (basically do what "Dissolve" in a GIS would do), but would like to avoid inserting "Merge Region" 60 times, opening the input mask for each of them and change the class filter. How can I achieve this "Merge by Class" all at once?


texts:
I've made a plugin asking for login credentials to a DB. Once logged in, the plugin is loading some layers from the DB. In the DB there are both admin and "normal" users.

I've saved two styles to the DB from the UI as default style for their layer respectively (seems you have to have DB admin rights to be able to save them).

Once I've saved the styles to the db I'm not able to see where they are stored i.e., I can't see any new tables. (I expected the styles to be saved in a table)

Running QGIS and the plugin again (as admin user) the layer styles are loaded automatically.

But, running QGIS and the plugin again (as normal user rnd) they are not loaded automatically. I then try to load the styles from the DB (from the GUI), but no styles are loaded.

Two questions.

How are the layer styles stored in the DB?
Is there a way to load the layer styles without DB admin rights?


texts: Is there a way to link eVis with time manager in QGIS

texts:
I have a point vector layer created from way points of a GPS. I then used the eVis plugin to view geotagged photos for each of these point. I also used the time manager plugin to play through the point layer. Is there anyway that as the time manager shows each point that I can see the photo associated with that point?

Thanks for any suggestions!


texts: How do I split a line feature in segments of a set interval distance in Qgis?

texts:
In order to record some vegetation along a Dutch highway,
I want to create a line feature along said highway of about 25 kilometers.
The line feature must be split every 100 meters to create seperate features that can have individual values stored to them.

Instead of manually breaking up the line in segments, I hope there is a faster way to do this.

I am using Qgis 3.03 Girona


texts: Trying to make one file out of many Hillshade files

texts:
I am trying to make a hilshade file for all of Colorado but when I download the DEM data for the state there are over 25 files. I have converted them all to hillsahde, but now I want to work with just one file instead of 25 of them.

Does anyone know of a way to do this?


texts: Why are my raster mean and standard deviation -1?

texts:
I have a set of about 25 rasters whose descriptive statistics I have been examining. The data in the rasters were calculated in Python and the rasters were created using Python-GDAL.

In all but one raster, the mean, std. dev., minimum and maximum seem correct. In one raster, however, the minimum and maximum  values are correct (4.842 and 9.458 respectively) but the mean and standard deviation are -1 if I take a look in ArcCatalog and -1.#IND in QGIS.

I opened the raster in Python via GDAL and calculated the minimum and maximum using data[data>band.GetNoDataValue()].min() and data[data>band.GetNoDataValue()].max() respectively, and the values these functions returned were correct (as above). But the mean value via data[data>band.GetNoDataValue()].mean() was 2.905, which is clearly nonsense as this value does not occur within the minimum and maximum value limits for this raster. Any ideas as to what's going wrong?


texts:
For my undergrad GIS project I am looking at whether left turns only (on UK roads... ie. not crossing traffic) could be more efficient for a delivery driver in a 2.2km radius.

I know it probably won't be, but I'm doing it for the sake of it, to show that I can use the software to the assessors etc.

Using OS Open Roads as a network and ArcMap's network analyst extension, is there a straightforward way to assign different costs to left and right turns? (i.e. make right hand turns take as long as four left hand turns or similar).


texts: Google Earth Engine large reduce resolution and reproject fails to export to drive

texts:
I'm attempting to aggregate (reduceResolution) a higher resolution image (e.g. 30m global forest loss as in code below) to a lower resolution image (e.g. 1000m pixels), based on some criteria, with a result where pixels and projection match (reproject) that of another lower resolution image (e.g. MODIS product at 1000m pixels as in the code below).
The code below performs as I intended it to, except that I cannot export the output at a large scale (e.g. the "aoi" given below) without the following error:
Error: Reprojection output too large (10429x8294 pixels).
I can export successfully if I reduce the region size ("aoi"), and the result is exactly as I would like it to be, but at a small scale. However, I would like to perform this analysis at a global scale.
My question is: How do I perform a global scale reduceResolution from 30m to 1000m and export it to my Drive?
Note. I use other data for my actual analysis, which explains various masking steps and crsTransform. These data are used for reproducibility.
var aoi =
    ee.Geometry.MultiPolygon(
        [[[[32.66972656249995, -7.554997127090209],
           [32.66972656249995, -7.554997127090209],
           [32.66972656249995, -7.554997127090209],
           [32.66972656249995, -7.554997127090209]]],
         [[[27.440234374999957, -3.9702890004775973],
           [27.440234374999957, -7.903365066051611],
           [33.46074218749995, -7.903365066051611],
           [33.46074218749995, -3.9702890004775973]]]], null, false);
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A2').first())
    .select('EVI');

// import hansen data for up to 2018 (v1.6)
var gfl = ee.Image("UMD/hansen/global_forest_change_2018_v1_6")

// Get information about the modis projection.
var modisProjection = modis.projection();
print('modis projection:', modisProjection);

// Get the forest cover data at modis scale and projection.
var gfl = gfl.select("loss")
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(), //Unweighted reducers so pixels are included if their centroid is in the region and the image's mask is non-zero.
      maxPixels: 10000
    });

// convert all modis pixels where 0-90% of pixels are masked to 1, otherwise they become 2
var gfl_modis = gfl.updateMask(gfl.mask().gte(0.1)).unmask(2)

Map.addLayer(gfl_modis, {min: 1, max: 2,  palette: ['DFDFDF', '171717']}, 'gfl_modis');

//Export the image, specifying scale and region.
Export.image.toDrive({
  image: gfl_modis,
  description: 'hansen_reproj',
  skipEmptyTiles: true,
  region: aoi,
  maxPixels: 1e13,
    shardSize: 1024,
  crs: modisProjection.wkt(),
  crsTransform: [1000, 0, -16449393.846139947, 0, -1000, 8745880.789507654]
});

Edit: I am aware of this question: Export of reprojected layer from GEE fails butbestEffort does not meet my requirements, and the answer doesn't allow me to make this analysis at a global scale.


texts:
I want to convert R 'library:raster' 'extent' object to a shapefile. I have the CRS available too.

  e=  c(4304916, 4305325, 365216, 365439)#xmin, xmax, ymin, ymax
  proj4string(e) = CRS("+proj=utm +zone=18 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

In older version of raster data there was a specific function ""
'extentTopolygon' which is depreciated now. My current approach is to create a raster object change its projection and assign the current extent and convert it to shapefile, which is not very efficient.


texts: Look Up table from (Latitude,Longitude) to country, city etc. boundaries (Bounding Box)

texts:
I'm looking for "Standard" assume file (?) containing best practice used simple "approximately squared" boundaries of countries, major cities etc.

The use case is getting (Lat,Lon) and perform look up search versus this static file.

If there is such ESRI shapefile etc. to put in PostGIS DB and then perform some query for this - That's also valid.


texts: Erasing boundaries between polygons that are next to each other using QGIS?

texts:
I have files that representes individual houses and blocks. I want to have a file that represents a residential area in a block.

Is there a way to erase the line, boundaries, between polygons that are adjacent?


texts: OpenLayers vector layer remove all features

texts:
I have a map (using OpenLayers) and a vector layer containing features(points). Sometimes I need to remove all features from the vector layer. When I use

vectorLayer.removeAllFeatures()

then it removes all features from only one zoom level. When I change zoom level, I see all my features again. How can I remove all features from all zoom levels?


texts: Looping dates in a feature class comparing them to date 1 month from now in ArcPy?

texts:
I thought this would be quite a simple task but I'm struggling. I want to loop over a date field in feature class looking for dates 1 month from today. It appears that I'm running into a simple date format issue but not sure what specifics to ask Google to figure it out for myself.

    #Import Modules
    import arcpy
    import datetime
    from dateutil.relativedelta import relativedelta

    #Today's date/one month from now
    today = datetime.datetime.now()
    checkDate = today + relativedelta(months=+1)

    #Featureclass variable
    fc = "C:\\Users\\99_Misc.gdb\\Check_v0"

    #Loop over expiry dates subtracting checkDate from date
    with arcpy.da.SearchCursor(extents,'Expiry_Date') as cursor:
        for row in cursor:
            row - checkDate

Traceback (most recent call last):
  File "C:/Users/Desktop/test.py", line 24, in <module>
    row - checkDate
TypeError: unsupported operand type(s) for -: 'tuple' and 'datetime.datetime'

Process finished with exit code 1

Do I need to make my date types consistent in some way?


texts: Benefit of DGPS for relative accuracy and speed measurements for receivers close-by each other

texts:
Consider that we are interested only in the speed and relative position accuracy of GPS receivers (i.e. their absolute try position on earth is NOT important). These receivers are all in a  90m x 90m area.
In this case, does it really matter whether we use DGPS or not?
For relative position/location accuracy, I think DGPS won't help much because the errors it corrects affect the receivers in a similar fashion (assuming that they see the same satellites, which is a reasonable assumption?)
Does DGPS help improve speed accuracy (assuming Doppler based speed)?


texts:
How can I get total feature count  of vector layer:

var mylayer = new OpenLayers.Layer.Vector(null, {
        "protocol": myWFSProtocol,
        "styleMap": style,
        "strategies": [
            new OpenLayers.Strategy.BBOX(),
            new OpenLayers.Strategy.Cluster()
        ]
        "filter": new OpenLayers.Filter.Logical({
            type: OpenLayers.Filter.Logical.OR,
            filters: [
                new OpenLayers.Filter.Comparison({
                    "type": OpenLayers.Filter.Comparison.NOT_EQUAL_TO,
                    "property": attributeNames.catid,
                    "value": "null"
                })
            ]
        })
    });

when I get mylayer.features.length it returns 5 but it actually 12. Is there any property of vector layer that gives total count.


texts: What are best practices to compare densities and locations between two sets of points?

texts:
I have two sets of point shapefiles.  They represent two different customer groups for a retail store.  Shapefile A represents online customers, Shapefile B represents in-store customers.

What I would like to do is create some sort of heat map that displays where the concentrations of these customers are.  That way I can hopefully show the differences in where these two different sets of customers are located in comparison to our actual store locations.  I even have how much money each individual customer has spent attributed in both shapefiles if I wanted to organize this by amount spent.  What do you think would be the best methodology to show this?  I have the Spatial Analyst and Business Analyst extensions installed.


texts: Keeping DEV, UAT, and PROD environments in sync with ArcSDE for Oracle?

texts:
Recently, I was doing some development work against our DEV environment. It seems logical at the time. When nearly complete, and after much face palming, I found the DEV environment was way out of sync with PROD (and UAT). Field names were different, nullable was true in PROD but not dev, data types were different, etc. etc.

What tools and practices could / should we be using to keep our various environments synchronized?

For example, in a pure Microsoft .NET & SQL Server world, I would use Red Gate to sync the database structure and data between the three environments.

What are my options for keeping ESRI SDE / Oracle synchronized between environments?


texts: arcgis server 10, vb.net, IMathOp problem

texts:
i am using arcgis server 10 and visual studio with vb.net to build web mapping application. inside my code i have the following line:

Dim RasterCal As ESRI.ArcGIS.SpatialAnalyst.IMathOp = _
   New ESRI.ArcGIS.SpatialAnalyst.RasterMathOpsClass

or even if i do it this way

Dim RasterCal As ESRI.ArcGIS.SpatialAnalyst.IMathOp = _
   New ESRI.ArcGIS.SpatialAnalyst.RasterMapAlgebraOp

in both cases i get initialization error which says:

The runtime application type must be specified before license
 initialization.

plz can you tell me what is wrong with the above examples of lines? thanks!


texts: Is there a way to automatically increment the Version value in the Config.esriaddinx file?

texts: I don't think there is an out-of-the-esri-box way to do this, but has anyone hacked up a way to change this value automatically for a build? Perhaps using a pre-build event command?


texts: How to remove warning message that appears when the src is not defined with pyQgis?

texts: I am currently learning to use PyQGIS, when loading a shapefile I get a warning message "The SRC is not defined", but I do not know what to do to prevent that message from appearing.


texts: Is age based SLD styling possible in GeoServer?

texts:
I have a set of features in a PostGIS table, where each feature has an attribute describing when the feature was created.

I would like to serve the features as a WMS layer and I would like to style each feature according to how old the feature is. An example could be a feature created today should be RED, if the feature is younger than 3 days old it should be GREEN and features older than that should be BLUE.

I could create a PostGIS view calculating the age of the features and style according to that, but it will also appear in the GetFeatureInfo lookup, which I would like to avoid.

Can this be done currently in GeoServer?


texts: Converting UTM32-EUREF89 shapefiles to align with Google maps

texts:
Im from Denmark, where the maps we can get from official hand is in UTM32-EUREF89 shapefile format. When these are applied on top of either google/oracle/bing they seem to be like 20-25 km of towards the south or the y-axis if you prefer.
Im using QGIS 2.6.1 and have tried converting the maps to 3857 and 4326 but that does nothing. I seem to be doing something wrong. Can anybody maybe give me a hint, as to where I should look.

Regards, Jacob


texts:
How can I build a feature request expression using variable inside expression?
For example

values = [1, 2, 3]
for j in range( y ):
it = vl.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"field_name" = values[ j ]' ) )

or (I have also tried sthing like this)

it = vl.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"field_name" = values[%d]' % j) )

Where could someone study how to build these kind of expressions?


texts: SCP QGIS Sentinel 2 plugin authorisation error

texts: I receive an http 401 auth failed error when trying to download Sentinel 2 images using the SCP plugin in QGIS. I registered a new account and confirmed my email but there is no support for new accounts after 21/12/2015. It seems to be possible to use a guest account but how does this work? . I tried this tutorial


texts: How to calculate the distance from the point to the raster' boundary in qgis

texts:
As the picture 1 shows,I use the points generate the Voronoi diagrams and change the polygons to the raster.The purpose of me is to calculate the min,max and mean distance from points to the raster'boundary they above(like picture 2).So that I can calculate the variance of every unit.

I do this in order to know the reachability of each grid in the area to the points.The greater the distance and variance, the worse the reachability.Can you help me? Or you have other ways?


texts: Using Python Opensource Libraries to check users who visited a store

texts:
I have a few files with location information:

1) User Lat Long Pings in Parquet Files (can be converted to CSV) format - this contains userid, timestamp, lat and lon
2) US County/State Shapefile
3) Store Centroid Data - centroid of stores like Walmart

Our goal is to:
1) Add geo attributes like county and state to user pings from county shapefile
2) Identify users who visited the store

Challenges:
1) Scale - User data has 30 million rows per day (over 6 months)

Can you please guide me on which geospatial libraries to use in order to create an efficient workflow? I briefly explored Fiona, Shapely, Geopandas.

I need to import the CSVs, shapefiles; define projections; create buffers; intersect them to get attributes and then output them. And I need to be able to do this at scale

Please let me know if my question is too vague


texts: Blur all other area outside polygon in mapbox-gl

texts:
I am trying to create an effect like

Can view it by adding feature from fiddle in
http://geojson.io/#map=14/46.9597/7.1446

when added the feature in mapbox as geojson layer of mapbox-gl I am getting the following output:

Fiddle for reference:https://jsfiddle.net/sanketg/24vzxcgo/

I am aware that according to mapbox v2 specification holes may cause problem in rendering pbf.

I have no idea regarding Geojson features.

What am I missing here?Is it that mapbox is not able to render the GeoJson feature neatly?

I am appending world bbox to polygon of GeoJson feature,to create the effect.Any other suggestions for better implementation are also welcome?


texts:
I'm trying to convert a bunch of shapefiles into geoJSON and found what looks like an easy way to do it but I can't get the code working. It doesn't matter if I run the code in OSGeo4W or just via idle I get the same error. (I do have OSGeo4W added to the system path.) My code and error are below. Could I get some guidance on what might be causing this?

import gdaltools

layer = <myFilePath>

ogr = gdaltools.ogr2ogr()
ogr.set_encoding("UTF-8")
ogr.set_input(layer, srs="EPSG:4326")
ogr.set_output('layer.geojson')
ogr.execute()

And the error:

Traceback (most recent call last):
  File "C:\Users\JOC-001\Desktop\Python Scripts\Random\Folium\GeoJSON\Con2GeoJSON.py", line 9, in <module>
    ogr.execute()
  File "C:\Python27\ArcGIS10.2\lib\site-packages\gdaltools\ogr2ogrcmd.py", line 289, in execute
    return self._do_execute(args)
  File "C:\Python27\ArcGIS10.2\lib\site-packages\gdaltools\basetypes.py", line 90, in _do_execute
    p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, bufsize=-1)
  File "C:\Python27\ArcGIS10.2\lib\subprocess.py", line 711, in __init__
    errread, errwrite)
  File "C:\Python27\ArcGIS10.2\lib\subprocess.py", line 948, in _execute_child
    startupinfo)
WindowsError: [Error 2] The system cannot find the file specified


texts: Determining coordinates of corners of raster layer using PyQGIS?

texts:
How can we determine the coordinates of the corners of the raster layer?

The image opens function QgsRasterLayer.


texts: Raster files are not exact in extent after clipping?

texts:
My two layers (chlorophyll-a and Bathmetry) are not exact after I clip my rasters. My bathymetry data (1) is reporting up to 10 decimal places where as chlorophyll-a (2) is only 3 decimal places and there is a rounding error occurring (see images). I need these two layers to be absolutely the same in extent, as i am uploading them as ascii files to maxent (species distribution), which requires exact extent. the extent i need is 118W, -10N, -18S and 126E.

I have tried clipping files by the extent of the other and also tried manually entering the values but after that the rounding error results in slight differences between the two files.

They are the same co-ordiante system and are the same cell size.

I am out of ideas to fix this can anyone suggest anything to solve this problem

The program is ArcMap 10.1


texts: download vector tile map from mapbox studio with python

texts:
I have a rendered mapbox map which is made up entirely by a few vector tileset. I wish to capture snipperts or AOI of the rendered map defined by users(or me) at varied zoom level into aa report
How should I go about doing it?
First way I have thought of is by trying to download the vector tiles using mapbox vector api. But it seem to be not right as I would need to download a few vector tileset to get lines and polygon (and still not the rendered map image!)
Other way I have thought of is through requesting the URL defined by WMTS/WMS service. Am I on the right track and how and where I could read up on help?
Thanks


texts: Merging two or more sub-basins into one?

texts:
I've been using Arc Hydro to analyze one of the unknown hilly regions with many watersheds. Unfortunately, Arc Hydro generates so many sub-basins that it looks really unrealistic and weird. However, I know that basins can be merged using the 'Basin processing toolbox'.

Does anyone know based on what criteria contiguous sub-basins are merged?


texts: Modify streaming interval in Collector for ArcGIS Online

texts: I'm using ESRI Collector App. I'm tracking location. Is there a way to modify the ESRI location tracking layer for it to use an interval under 5 seconds?


texts: Making vertical spacing between parameters smaller on ArcGIS Script tool dialog

texts: I would like to make space between parameters smaller in ArcGis Script tool. For example, in the image attached I would like to make the distance between two parameter Standard terrain and User defined smaller. These parameters are Boolean data type.


texts: Can't create named map containing complex CSS condition using cURL

texts:
I'd like to style contour lines so that they show differently at varying levels of zoom.  When zoomed right out, only the 200m, 400m, 600m, etc. lines should show.  When zoomed in a bit further, the 100m, 300m, etc. lines should appear.  Text should appear too when zooming in.

I can achieve this fine using reasonably concise OR syntax in Builder, but when I try to cURL definitions containing OR conditions to create a named map, the map stops working.

For example, this one works both when cURLed and in Builder:

#layer[value='200'] {
    line-width: 0.7;
    line-color: #4b3ebc;
    line-opacity: 1;
}

But this one containing the OR condition ,[value='400'] only works in Builder: it does not work when cURLed:

#layer[value='200'],[value='400'] {
    line-width: 0.7;
    line-color: #4b3ebc;
    line-opacity: 1;
}

I can achieve the desired effect by cURLing this more verbose expression of the same thing, but would rather use more concise syntax because the full styling I need involves many more conditions:

#layer[value='200'] {
  line-width: 1.5;
  line-color: #4b3ebc;
  line-opacity: 1;
}

#layer[value='400'] {
  line-width: 1.5;
  line-color: #4b3ebc;
  line-opacity: 1;
}

For reference, here's the full styling I'm trying to achieve, expressed using concise syntax:

#layer[value='200'],[value='400'],[value='600'],[value='800'],[value='1000'],[value='1200'],[value='1400'],[value='1600'],[value='1800'],[value='2000'] {
  line-width: 1.5;
  line-color: #4b3ebc;
  line-opacity: 1;
}

#layer[zoom>10][value='200'],[zoom>10][value='400'],[zoom>10][value='600'],[zoom>10][value='800'],[zoom>10][value='1000'],[zoom>10][value='1200'],[zoom>10][value='1400'],[zoom>10][value='1600'],[zoom>10][value='1800'],[zoom>10][value='2000'] {
  text-face-name: 'DejaVu Sans Book';
  text-name: [value];
}

#layer[zoom>10][value='200'],[zoom>10][value='400'],[zoom>10][value='600'],[zoom>10][value='800'],[zoom>10][value='1000'],[zoom>10][value='1200'],[zoom>10][value='1400'],[zoom>10][value='1600'],[zoom>10][value='1800'],[zoom>10][value='2000'] {
  text-face-name: 'DejaVu Sans Book';
  text-name: [value];
}

#layer[zoom>11][value='100'],[zoom>11][value='300'],[zoom>11][value='500'],[zoom>11][value='700'],[zoom>11][value='900'],[zoom>11][value='1100'],[zoom>11][value='1300'],[zoom>11][value='1500'],[zoom>11][value='1700'],[zoom>11][value='1900'],[zoom>11][value='2100'] {
  line-width: 1.0;
  line-color: #4b3ebc;
  line-opacity: 1;
  text-face-name: 'DejaVu Sans Book';
  text-name: [value];
}

#layer[zoom>12][value='50'],[zoom>12][value='150'],[zoom>12][value='250'],[zoom>12][value='350'],[zoom>12][value='450'],[zoom>12][value='550'],[zoom>12][value='650'],[zoom>12][value='750'],[zoom>12][value='850'],[zoom>12][value='950'],[zoom>12][value='1050'],[zoom>12][value='1150'],[zoom>12][value='1250'],[zoom>12][value='1350'],[zoom>12][value='1450'],[zoom>12][value='1550'],[zoom>12][value='1650'],[zoom>12][value='1750'],[zoom>12][value='1850'],[zoom>12][value='1950'],[zoom>12][value='2050'],[zoom>12][value='2150'] {
  line-width: 0.7;
  line-color: #4b3ebc;
  line-opacity: 1;
  text-face-name: 'DejaVu Sans Book';
  text-name: [value];
}

I won't post the same CSS expressed the verbose way - if you consider the way those OR conditions decompose, clearly it's a lot more repetitive.

(I should note that in fact I use Postman, not cURL, but all the docs talk about cURL, so I thought I'd use the same language here.)


texts: GetFeatureInfo Request with GeoServer and OpenLayers.loadURL Not Working

texts:
I'm making a simple GetFeatureInfo request against a GeoServer WMS, but it returns empty responses. The request utility I'm using is OpenLayers.loadURL.

The funny thing is, I can paste the generated request URL from FireBug into a new browser tab and get a valid response. So, it seems to be something with the OpenLayers.loadURL. See the examples below.

I can get a valid response with this map event where the GetFeatureInfo request is in a new window:

  map.events.register('click', map, function (e) {
    var url = "http://localhost:8080/geoserver/MYLAYER/wms"
    var params = ""
      + "REQUEST=GetFeatureInfo"
      + "&SERVICE=WMS"
      + "&EXCEPTIONS=application/vnd.ogc.se_xml"
      + "&VERSION=1.1.1"
      + "&BBOX=" + map.getExtent().toBBOX()
      + "&X=" + Math.round(e.xy.x)
      + "&Y=" + Math.round(e.xy.y)
      + "&INFO_FORMAT=text/html"
      + "&QUERY_LAYERS=NERD:NERD Risk Rates"
      + "&LAYERS=NERD:NERD Risk Rates"
      + "&FEATURE_COUNT=50"
      + "&format=image/png"
      + "&SRS=EPSG:900913"
      + "&STYLES="
      + "&WIDTH=" + map.size.w
      + "&HEIGHT=" + map.size.h;

    console.log(url+"?"+params);
    window.open(url+"?"+params,
      "getfeatureinfo",
      "location=0,status=0,scrollbars=1,width=600,height=150"
    );

   });

However, this map event using OpenLayers.loadURL always returns an empty response and no error. But if I cut and paste the GET request url from FireBug into a new browser window, I get the proper response:

  map.events.register('click', map, function (e) {
    var url = "http://localhost:8080/geoserver/MYLAYER/wms"

    var paramHash =  {
      REQUEST : "GetFeatureInfo" ,
      EXCEPTIONS : "application/vnd.ogc.se_xml" ,
      SERVICE : "WMS" ,
      VERSION : "1.1.1" ,
      BBOX : map.getExtent().toBBOX() ,
      X : Math.round(e.xy.x) ,
      Y : Math.round(e.xy.y) ,
      INFO_FORMAT : "text/html" ,
      QUERY_LAYERS : "NERD:NERD Risk Rates" ,
      LAYERS : "NERD:NERD Risk Rates" ,
      FEATURE_COUNT : 50 ,
      format : "image/png" ,
      srs : "EPSG:900913" ,
      styles : "" ,
      WIDTH : map.size.w ,
      HEIGHT : map.size.h
    };

    OpenLayers.loadURL(url, paramHash, this, setHTML, setHTML);
        OpenLayers.Event.stop(e);

   });

Someone have any ideas?


texts: Determining Personal Geodatabase version?

texts:
I am using ArcMap software.

How can I find out Personal Geodatabase Version?

I have a client file but I want to find out which version they have used.


texts: While measuring an angle in QGIS using measure angle tool, it is showing as 0 to 180 degree for clockwise direction and -1 to -179 for anticlockwise. What I need is, interior angle between vertices of irregular polygon. How should I calculate in QGIS?


texts: Google Earth Engine, 'computation time out' error when calculating area of binary masked

texts:
Using Landsat 8 TOA, I tried to find the area of mndwi pixel using binary masking, but its giving computation time out.
https://code.earthengine.google.co.in/c1a91260ff1628376dd19ff7f938261d
asset:https://code.earthengine.google.co.in/?asset=users/013limmazumms/Maharashtra_district


texts:
I am trying to convert a geometry column to a geography column and calculate the area using this query:

SELECT dccode,SUM(CAST(St_Area(geom::geography,true) AS NUMERIC)) as Area
       FROM  schema.table
       GROUP BY dccode

but I am getting a strange error:

...
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY

ERROR:  lwgeom_area_spher(oid) returned area < 0.0
********** Error **********

ERROR: lwgeom_area_spher(oid) returned area < 0.0

I've run this query to verify that there is no invalid geometry:

select gid, dccode, st_isvalidreason(geom) from table.schema where st_isvalid(geom) = false;

but it doesn't return anything. How can I isolate the geometry that is < 0.0?


texts: Georeferencing files that are orthogonal to a point on a sphere

texts:
I have about 65000 photos, each of them a square where an inscribed circle is radius 1 degree.

Each of these photos has an associated latitude and longitude of the center.

I'm trying to figure out how to georeference these photos with gdal so that I can warp them into a large, more sane latlong coordinate system.  Once I figure out one of them, it will be trivial to automate the rest, but I just can't figure out how to properly georeference an orthogonal photo.

For instance, tile 1 has a central latitude and longitude of 0.  Further, along the central X-axis it occupies 1 deg and along the y-axis, it occupies 1 degree.  Using the attached photo as a guide:

NOTE: (long, lat)
A: 0,0
B: 0,.5
C: .5,.5
D: 0,-.5
E:-.5,-.5

The next photo in this series along longitude would have the following coordinates:

A:1,0
B: 1,.5
C: 1.5,.5
D: 1,-.5
E:.5,-.5

--

Likewise, a photo at the pole would have a central latitude of 90 and longitude of 0.  So

A: 0,90
B: 0,89.5
C: 90,89.5
D: 180,89.5
E: -90,89.5

Further, the next photo in this series would be rotated by 1 degree.  So

A: 1,90
B: 1,89.5
C: 91,89.5
D: -179,89.5
E: -89,89.5

An so and, so forth.

--

My proposed solution is something like:

gdal_translate -a_srs "+proj=ortho +lat_0=0 +lon_0=0 +a=6378140 +b=6378140 +units=d" input_000_00.jpg output_000_00.tif

But, when I do that and attempt to warp it, I get the following:

  ERROR 1: Unable to compute a transformation between pixel/line and georeferenced coordinates for test17990.tif.  There is no affine transformation and no GCPs.


texts: Is the Geometry(4326) an alias for Geography datatype?

texts:
The main reference about Geography datatype seems docs/using_postgis_dbmanagement:

Prior to PostGIS 2.2, the geography type only supported WGS 84 long lat (SRID:4326). For PostGIS 2.2 and above, any long/lat based spatial reference system defined in the spatial_ref_sys table can be used (...)
Regardless which spatial reference system you use, the units returned by the measurement (...)  are in meters.

But it is not explicit: when I am "casting to another type" and when PostGIS is using internally the same "Geography calculations", without any casting need?
Other texts, like postgis-intro/geography, and other documents are also confuse: ... is a Geometry with no projection a Geography? The "standard no projection" is Geometry(4326), so, why the need for a Geography data type if PostGIS can detect it by SRID?
Of course, for PostgreSQL is important to change the data type, to do function overloading... But the guide not explains the equivalence schemas in the internal representation.

NOTES
There are some previous answers to earlier questions:

this answer to Oracle equivalent to postgis EPSG 4326 "geometry" type, and
this answer to Understanding geography, geometry, distance and wkt in PostGIS

but those questions are not duplicates of this one.
About comments, to avoid terminological confusion: a Geometry with SRID=4326 is not "only lat long". The so-called "WGS84" is a CRS (coordinate reference system), a synonym  and human-readable abbreviation for urn:ogc:def:crs:EPSG::4326. Any formal standard CRS is composed by   ellipsoid parametrized model reference + parametrized datum reference, so, the formal definition is:

the WGS84's ellipsoid  is urn:ogc:def:ellipsoid:EPSG::7030;
the WGS84's  datum is urn:ogc:def:datum:EPSG::6326.


texts: Alter Fields in multiple feature classes in multiple workspaces

texts:
I have 22 geodatabases all in one folder. Each gdb has two feature classes. For each feature class, I need to rename a few of the fields. I've had success with this code, but it only works for each GDB at a time:

import arcpy
from arcpy import env
env.workspace = r'C:\...\filename.GDB'

for fc in arcpy.ListFeatureClasses('*'):
    try:
        arcpy.AlterField_management(fc, 'ADDR_SN', 'STREET_NAME', 'STREET_NAME')
    except:
        pass

As I said, that code is just for one GDB at a time, so I'd have to change the name in env.workspace each time. That's not a big deal, I just want to see if it's possible to iterate through every GDB in the folder. So I tried adding into the code:

for workspace in arcpy.ListWorkspaces('*', 'All'):

When I add the above code, I also change the env.workspace to the folder name that contains all the geodatabases, but nothing happens. No error is returned, nothing happens.


texts: QGIS 3: About "topological editing" and 'avoid intersection"

texts: Where are those two editing options ? They used to be with the snapping options but apparently they have been moved somewhere else...


texts: Scaling two rasters to same elevation scale in QGIS?

texts:
I'm pretty new to using GIS, but I've got two overlapping rasters of the same general area - one's a 1m resolution DEM from the Pennsylvania side, and the other's a 0.6m DEM resolution from the Maryland side. Both have their elevations measured in feet, but the elevations from the two DEMs in the overlapping areas are drastically different. Here's a screenshot from trying to display them on the same elevation scale: https://i.imgur.com/PgoNJwK.png

Here's the spatial reference information for the PA side:

Geodetic_Model:

    Horizontal_Datum_Name:North American Datum of 1983
    Ellipsoid_Name:Geodetic Reference System 80
    Semi-major_Axis:6378137.000000
    Denominator_of_Flattening_Ratio:298.257222

Vertical_Coordinate_System_Definition:

Altitude_System_Definition:

    Altitude_Datum_Name:North American Vertical Datum of 1988
    Altitude_Resolution:0.00000
    Altitude_Distance_Units:feet
    Altitude_Encoding_Method:Explicit elevation coordinate included with horizontal coordinates

And the corresponding CRS in QGIS for the PA side: +proj=lcc +lat_1=39.93333333333333 +lat_2=40.96666666666667 +lat_0=39.33333333333334 +lon_0=-77.75 +x_0=600000 +y_0=0 +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

And the MD side:

Geodetic_Model:

    Horizontal_Datum_Name:North American Datum of 1983
    Ellipsoid_Name:Geodetic Reference System 80
    Semi-major_Axis:6378137.000000
    Denominator_of_Flattening_Ratio:298.257

Vertical_Coordinate_System_Definition:

Altitude_System_Definition:

    Altitude_Datum_Name:North American Vertical Datum of 1988 (GEOID12A)
    Altitude_Resolution:0.0
    Altitude_Distance_Units:feet
    Altitude_Encoding_Method:Implicit coordinate

And the MD CRS: +proj=lcc +lat_1=39.45 +lat_2=38.3 +lat_0=37.66666666666666 +lon_0=-77 +x_0=399999.9999999999 +y_0=0 +datum=NAD83 +units=us-ft +no_defs


texts: Creating Multipart Line from coordinate pair lists using Geometry objects in ArcPy?

texts:
I am developing a script in which I want to create a multipart line geometry.

For example:

I have one array of points a,b,c and second array of points e,d,f.

I want to create one multipart line where a,b,c defines one part and e,d,f defines the second part.

Is this possible using ArcPy?

Something like arcpy.polyline([a,b,c], [d,e,f]) ?


texts: ogr2ogr window closes after termination?

texts:
I'm using ogr2ogr for export from PostGIS to FileGDB. It works fine, but after termination the command window closes immediately.

Is there a way to keep it open and get some kind of message or report if the task was executed successfully?


texts:
I'm trying to create a buffer around a multi-linestring object with st_buffer.

Here's my initial attempt:

library(rnaturalearth)
library(dplyr)
library(sf)
library(rmapshaper)

proj_sahel <- "+proj=aea +lat_1=9.71 +lat_2=22.58 +lon_0=2.81"

boundaries <- ne_download(scale = 10, category = "cultural", type = "admin_0_boundary_lines_land", returnclass = "sf")

buffer50 <- boundaries %>% filter(adm0_a3_l == "MLI" & adm0_a3_r %in% c("BFA", "NER")) %>%
st_transform(crs = proj_sahel) %>%
st_buffer(dist = 50000) %>%
st_cast("POLYGON")

It works fine except for the unfortunate artifact right in the middle of the buffer zone as you can see below:

I tried simplifying the geometry before and after st_transform with ms_simplify, it didn't make any difference. How can I avoid this or clean it after the fact?


texts:
Problem:

I am using ogr2ogr to convert osm file to shp and some features realated to the lakes are missing. I noticed for ex. lake is missing when it is flow-through lake (river is at both ends).

My result

Shapefile with some lakes but many are missing.

Question

Can someone tell me what I can do to get those missing lakes?

Part of my Osmconfig.ini

closed_ways_are_polygons=natural=water

[multipolygons]

attributes=natural

Picture:

You can only see river polygon here but there was no point to add any more layers, because I am interested in river and lakes.

rose polygon = river polygon

basemap = OSM Basemap in ArcGIS Pro

Edit:

I found out that when lake is within 'relation' it is correctly converted to the record in the shp, but when lake is within 'way' it is missing in shp.

Version:

GDAL 2.4.0, released 2018/12/14


texts: How to downsample raster image for nightlights data in Python?

texts:
I'm working with the nightlights data - 30 arc second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude.
I want to downsample it to 5 arc-min resolution.

Now my question is how to best do it? I'm using rasterio to load the data and looked at this way but it does not work for me as I'm using python 3 and get error - "TypeError: 'float' object cannot be interpreted as an integer" I'm also worried as to how it will impact my coordinate systems etc.

I haven't found anything latest on this yet. What's the best way to go about it?

nlight_path = "./data/nightlights/F182013.v4c_web.stable_lights.avg_vis.tif"
# open raster data
with rio.open(nlight_path) as src:
    nlight = src.read(masked=True, indexes = 1)
    nlight_meta = src.profile
    nlight_extent = rio.plot.plotting_extent(src)
print(nlight_extent)
print(nlight_meta)
print(nlight.shape)
print(type(nlight))

(-180.00416666665, 180.00416522665, -65.00416610665, 75.00416666665)
{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 43201, 'height': 16801, 'count': 1, 'crs': CRS.from_dict(init='epsg:4326'), 'transform': Affine(0.0083333333, 0.0, -180.00416666665,
       0.0, -0.0083333333, 75.00416666665), 'tiled': False, 'interleave': 'band'}
(16801, 43201)
<class 'numpy.ma.core.MaskedArray'>

Also can I directly do the downsampling on masked array? But then how will it affect the mask?


texts: Rasterize solid 3D-building blocks with QGIS

texts:
Stockholm city provides 3D building blocks as open vector data. I want to rasterize these solid cubes to be able to perform a multi criteria analysis. But no matter which rasterize-command I try in QGIS, the output .tiff-file turns black!? The Z-column in my vector file is called "Z".

Any ideas?

This is how the attribute data for the 3D blocks looks like:


texts: TopoJSON coordinates are off with actual coordinate using Leaflet

texts:
I have this TopoJSON file that I have been working on quite sometime now. and recently I tried plotting it using Leaflet but the polygons are off with the actual area it surrounds. Here's the image.

Is there a way to adjust the polygons with their actual position within the map? I have a shapefile of it, can I do it with QGIS? Or there's a third-party app that I can use to adjust it?


texts: Printing/Exporting maps based on contour lines using ArcGIS Pro?

texts:
I have generated a map with layers for contour lines, roads, trails, rivers, lakes, etc. What I would like to do is print to PDF (AI or other) a page/map for each contour level. Ex: one for 1900m with the features that lie on that contour, one for 2000m with those features and so forth. Think of it like a flip book where each page goes up 100m and shows the associated features.
Here is an example of the first and second contours and the information present on them.

Currently I am printing to SVG format and then manually separating everything in Adobe Illustrator... but there has to be a better way!

I have access to ArcGIS Pro, although am more skilled in QGIS.


texts:
I'm working with Sentinel-2 satellite and i'm working with 2,3,4,8 bands. I have to make composite mergining bands. I download images from https://scihub.copernicus.eu/. I used gdal_merge in console for creating composites. I have batch files, code below:

For /F  %%i in ('where /R %cd% *B02*.jp2') do set file1=%%i
 For /F  %%n in ('where /R %cd% *B03*.jp2') do set file2=%%n
 For /F  %%m in ('where /R %cd% *B04*.jp2') do set file3=%%m
 For /F  %%o in ('where /R %cd% *B08*.jp2') do set file4=%%o
gdal_merge -o composite.tif %file1% %file2% %file3% %file4% -co COMPRESS=LZW -separate
pause

It's working, but i'd like to write own script using Python.

Could somebody give me advice how can i do it, please ? Is there another way to do it ?


texts: Better way to handle HDF and GeoTIFF files

texts:
I am using GDAL and MODIS reprojection tool in my python project and want to know better way for math new indices of HDF data (MODIS data).
Also I want to visualize them in web. Which format is more flexible for math and visualize my results?


texts: Add Undo/Redo functionality in ArcMap to ArcPy script with manually started edit session

texts:
I've written a set of simple ArcPy scripts for use with ArcMap 10.5.1 in versioned PostgreSQL Enterprise geodatabases. They're loaded into an add-in toolbar and intended to be used like a speed dial of shortcuts for updating specific attributes during a manually started edit session within ArcMap.
import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
SelectedLayer = [lyr.longName
    for lyr
    in arcpy.mapping.ListLayers(mxd)
    if lyr.getSelectionSet()]
#Stops script unless there is a selection
if len(SelectedLayer) < 1:
    print "Selection required"
    exit()
#Iterates through selected layers and applies 1 to all selected features for Test attribute
for layer in SelectedLayer:
    ScaleLayer = arcpy.mapping.Layer(layer)
    ScaleSelection = ScaleLayer.getSelectionSet()
    with arcpy.da.UpdateCursor(ScaleLayer, 'Test') as ucursor:
        for row in ucursor:
            row[0] = 1
            ucursor.updateRow(row)
arcpy.RefreshActiveView()

The script works, but I've hit a roadblock for including Undo/Redo functionality, which is vital.
startOperation() and stopOperation() seem to be the way to do this, but the commands don't recognize that the user is already in an edit session. To use these commands, I have to also add startEditing, which breaks the Undo/Redo stack in place and can result in strange errors if Undo is used. Oddly, including stopEditing after stopOperation also does not recognize an edit session is in place. Here is a working version of the code that includes those elements:
import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
SelectedLayer = [lyr.longName
    for lyr
    in arcpy.mapping.ListLayers(mxd)
    if lyr.getSelectionSet()]
#Stops script unless there is a selection
if len(SelectedLayer) < 1:
    print "Selection required"
    exit()
#Iterates through selected layers and applies 1 to all selected features for Test attribute
for layer in SelectedLayer:
    ScaleLayer = arcpy.mapping.Layer(layer)
    workspace = ScaleLayer.workspacePath
    edit = arcpy.da.Editor(workspace)
    ScaleSelection = ScaleLayer.getSelectionSet()
    edit.startEditing(True, True)
    edit.startOperation()
    with arcpy.da.UpdateCursor(ScaleLayer, 'Test') as ucursor:
        for row in ucursor:
            row[0] = 1
            ucursor.updateRow(row)
    edit.stopOperation()
arcpy.RefreshActiveView()


texts: Avoid reshaped/jagged polygons when exporting from ArcMap to PDF

texts:
I realised that exporting a huge amount of building polygons (in my case 86,000 polygons) to a A2 PDF within ArcGIS 10.2.2 the polygons get unacceptably jagged or even reshaped. It gets a bit better if the dpi value is set to a higher number (e.g. 1000) but it's still not the desired quality.

Polygons within my ArcMap data frame:

Polygons after export to a PDF (400 dpi):

Is there a way to avoid this problem?


texts: Converting raster to polygon through ModelBuilder?

texts:
I have a set of raster which I want to convert to the polygon.
Below is the model I prepared, but its showing below error.

I checked the error code, which explains Alter the name to remove any invalid characters. For instance, consider using an underscore (_) in place of a space

I don't have any space in my input files and in the output file I am writing "%Name%".

What exactly I am doing wrong?

I am using ArcGIS 10.7.1 for this.


texts: QGIS label: function to get a value from an other layer

texts:
I have a layer called "Projet", which overlays a layer called "Parcels".

I have a symbol render which allows me to highlight the touched parcels.

The expression I use is:

Is it possible in in the label of my layer "Projet" to put the number of the parcels (attribut of the layer Parcels?

[EDIT]
So I wasn't clear enough...
As you can see in the first picture, each parcel has a unique NUMBER.
What I trying to do, is to get this unique number for each intersected parcels.
I try to have the number value for each centroid. On the example above (second picture), I have 4 parcels in red. So I would like 4 unique numbers, one for each centroids! Some thing like that:

Sorry for my bad english if it's not clear enough


texts: Edited features only store edits upon map refresh

texts:
I am setting editing toolbars for features on my map as follows:

 for (var graphicNum = 0; graphicNum < graphicList.length; graphicNum++) {
        if (graphicList[graphicNum].symbol != undefined) {
            //Initializing the toolbar to be used for the feature modification process
            moveToolbarLabel = new esri.toolbars.Edit(app.map);
            moveToolbarLabel.activate(esri.toolbars.Edit.MOVE | esri.toolbars.Edit.SCALE | esri.toolbars.Edit.ROTATE, graphicList[graphicNum]);
            allToolbars.push(moveToolbarLabel);
        }
    }

When moving or scaling a feature after the scale and rotation points appear the feature reverts to its original state when the mouse is unclicked (Meaning click on scale point, drag to resize, let go and feature goes back to original size). However if the map is then refreshed by zooming in the feature takes on the changes that would have been made when the mouse was dragged. Is there any way to make the feature keep the edits right away rather than having to zoom in/out of the page to see the changes?


texts:
I would like connect to ArcGIS server, but I have an error message
"la référence d'un objet n'est pas définie à une instance a une objet".
I have many examples:

string serverName = "adresse";
string instance = "ArcGIS/services";
// Using the ADF Common API data source implementation of the ArcGIS Server Web Service API
    string serviceNameAndType = "nameservice/MapServer";
    ESRI.ArcGIS.ADF.ArcGISServer.MapServerProxy agsSoapMapServerProxy = new
ESRI.ArcGIS.ADF.ArcGISServer.MapServerProxy("http://" + serverName + "/" + instance + "/" + serviceNameAndType);

or

 ESRI.ArcGIS.ADF.Identity identity = new ESRI.ArcGIS.ADF.Identity("usernamer", "password", "domain");
 ESRI.ArcGIS.ADF.Connection.AGS.AGSServerConnection agsconnection;
agsconnection = new ESRI.ArcGIS.ADF.Connection.AGS.AGSServerConnection("nomserver", identity);
 agsconnection.Connect();
 IServerObjectManager SOM = agsconnection.ServerObjectManager;

or

string mapServiceUrl = "http://servername/ArcGIS/services/servicename/MapServer";

ESRI.ArcGIS.ADF.ArcGISServer.MapServerProxy mapServerProxy = new ESRI.ArcGIS.ADF.ArcGISServer.MapServerProxy(mapServiceUrl);

I hope that someone can help me :)


texts:
I have a point where a Gpolyline line is drawn from a marker to another point. After I click on it a second time, I want to remove the polyline. How do I remove the line after clicking on the marker a second time?

My pseudo-code is as follows:

  google.maps.event.addListener(markersArray[0], 'click', function() {
     'draw line to point'
  } );

(I did also ask this question on StackOveflow - I'm not sure whether it is ok to post the same question here)


texts: Support viewing correctly and creating NSG GPKG RASTER TILES

texts:
Trying to view a US Navy Or NGA NSG Geopackage raster tiles in OpenLayers or Leafletjs webmaps
I have used NGA geopackage-js but it NOT been updated to support these NSG spec geopackage that do not use zoom level and column Row numbering like  Google XYZ TOP LEFT ORIGIN or TMS bottom left origin. In addition, they do nkt start with zoom level 0 being one tile. It uses two tiles. samples can be found at https://geoint.nrlssc.navy.mil/#/ or use their download area of interest tool to create a new geopackage. https://github.com/ngageoint/geopackage-js


texts: ArcGIS 10.1 using Repast's Agent Analyst Extension - refresh.exe does not work

texts: I have installed the agent analyst extension for ArcGis. The refresh.exe, however, does not work. I have made a direct reference to the folder where refresh.exe is.


texts: Display GeoPDF's coordinates in metres created using QGIS

texts:
I exported a GeoPDF map using QGIS 3.10.1. My project coordinates system is in BNG (EPSG 27700) and so the CRS in map grid properties in the layout. In other words, source data, the project CRS and the grid properties are in metres.

But the PDF I created earlier doesn't show the units in metres when it is opened using Adobe reader (Geospatial Location Tool). And I would like to see the coordinates of the GeoPDF in metres (Easting and Northing).


texts:
Here the code I have written to a shapefile in MapGuide Open Source (in PHP running on IIS 7.5). The code is actually the modified class MarkupEditor (markupeditor.php) from the General Tracks (GT) example.

class MarkupEditor
{
    private $args = null;
    private $site = null;
    function __construct($args){
        $this->args = $args;
        $this->site = new MgSiteConnection();
        $this->site->Open(new MgUserInformation($args['SESSION']));
    }
    //
    // Rest of the original code in the class
    //
    function AddPointToShp($featureID){
        $geometryFactory = new MgGeometryFactory();
        $agfWriter = new MgAgfReaderWriter();
        $vertices = explode(",", $this->args['GEOMETRY']);
        $coord = $geometryFactory->CreateCoordinateXY((double) $vertices[0], (double) $vertices[1]);
        $point = $geometryFactory->CreatePoint($coord);
        $point = $point->Transform($this->GetTransform());
        $byteReader = $agfWriter->Write($point);
        $propertyValues = new MgPropertyCollection();
        $propertyValues->Add(new MgInt32Property("FeatId", $featureID));
        $propertyValues->Add(new MgStringProperty("TYPE", "Breader"));
        $propertyValues->Add(new MgStringProperty("DESC", "OCB'));
        $propertyValues->Add(new MgGeometryProperty("Geometry", $byteReader));
        $this->InsertMarkupFeatureToShp($propertyValues);
    }

    function InsertMarkupFeatureToShp($propertyValues){
        try{
            $featureService = $this->site->CreateService(MgServiceType::FeatureService);
            $featureSourceId = $this->args['MARKUPLAYER'];
            $commands = new MgFeatureCommandCollection();
            $commands->Add(new MgInsertFeatures($this->args['OPENMAKUP'], $propertyValues));
            $featureService->UpdateFeatures($featureSourceId, $commands, false);
        }
        catch (MgFeatureServiceException $ex){
            echo '<b>MgFeatureServiceException : </b>'.$ex->GetExceptionMessage();
        }
        catch (MgInvalidArgumentException $ex){
            echo '<b>MgInvalidArgumentException : </b>'.$ex->GetExceptionMessage();
        }
        catch (MgInvalidOperationException $ex){
            echo '<b>MgInvalidOperationException : </b>'.$ex->GetExceptionMessage();
        }
        catch (MgFdoException $ex){
            echo '<b>MgFdoException : </b>'.$ex->GetExceptionMessage();
        }
    }
}

The members of $args array passed in the constructor are:

Array
(
    [MARKUPLAYER] => Library://DataStore/Data/Electrical/Breaker.FeatureSource
    [SESSION] => 970bab06-4bd2-11e1-8000-000000000000_en_7F0000010AFC0AFB0AFA
    [MAPNAME] => DataStore
    [OPENMARKUP] => Default:breaker
    [TGT] => 1
    [POPUP] => 0
    [GEOMETRY] => 1.5041909054501184,0.39827301781943014
)

$featureID in function AddPointToShp() is a numeric value. The code always stops execution at this line in function InsertMarkupFeatureToShp():

$featureService->UpdateFeatures($featureSourceId, $commands, false);

It does not even generate an exception so that I may be able to find out what is wrong.

Some posts on osgeo.org say that GT Example is used to write to SDF only. Where as, I am able to read any field from a desired shapefile using some other modified functions in the same file.

Any help please?

UPDATE

I have converted the target Coordinate System to the one as in existing shapefile. Now I get the MgInvalidArgumentException:

Invalid argument(s): [1] = " " The string cannot be empty.

Also I have tested $propertyValues [$propertyValues->GetCount()]. It returns the count as 4.


texts: How to make list of nodes at end points of polylines that connect them to others using ArcGIS for Desktop?

texts:
I have a sewer network represented as a connected pipelines in a GIS shapefile.

the network is very close to a tree graph that each node is only connected with one or two nodes.

I need to reconstruct the layout of this network on another program.

I first used vertices to points feature to specify the start and the end point of each line.

However, this doesn't help too much because I need to specify which node is connected to which node and the feature of the line that is connecting these two nodes like its diameter or its FID for example.


texts:
I have a table with 2000 records and with those columns
OutletName, longitude, latitude, geom

I would like to find a way to find the farthest outlet between a selection of outlets

For example we have outletnames : A,B,C,D,E,F,G,H,I,J,K...

And we would like to know who is the farthest outlet from A
between those outlets 'B' 'F' and 'D'

Now I am doing it with lot of queries calculating distances between all of them and comparing the distances in java, but it takes time. Is there a better way to do it?


texts:
I have a QGIS project file, it has a point layer in GPKG format and online map layers. Coordinate System of online map layers is WGS 84 / Pseudo-Mercator. Point layer's default coordinate system is WGS84 4326. I try to modify coordinate system of point layer using QGIS C++ API and QgsQuick based on user's input with EPSG code. The code below creates new project. When I open the project and create new points on online map layer, it is OK. But when I open the project again, the points are at wrong place, like middle of Africa. I think there is a problem about on-the-fly projection about the project file. I wrote this C++ code from an answer about creating project with PyQGIS.
Is there a way to fix this problem? Should I set the extent, project coordinate system before saving the project?
You can see the project's source code here.
QString ProjectModel::addNewProject(QString name, QString epsgcode )
{
    // create CRS with EPSG
    QgsCoordinateReferenceSystem my_crs = QgsCoordinateReferenceSystem( "EPSG:" + epsgcode );
    if ( !(my_crs.isValid()) ) {
        return "crs_not_found";
    }

    // new project folder path name
    QString new_path = dataDir() + "/" + name;

    // folder path name to directory
    QDir dir( new_path );

    if ( !(dir.exists()) ) {
        // Copy sample qgs project and point layer to destination path
        copy_survey_project( name );
        // New project path
        QString projectName = new_path + "/" + name + ".qgs";
        // initialize project and clear project
        mProject = QgsProject::instance();
        mProject->clear();
        // Read the copied qgs project
        mProject->read( projectName );
        // Point layer
        QgsMapLayer *myLayer = mProject->mapLayersByName("survey")[0];

        if( myLayer && myLayer->isValid() ){
            // Set new CRS to the layer
            myLayer->setCrs( my_crs );
            // Save the project
            mProject->write( projectName );
            mProject->clear();
        }
        return "ok";
    }
    else{
        return "fileExists";
    }
}


texts: Why does a raster that consists of 3 bands appear greyed out?

texts: I have a raster that is formed by stacking the bands blue, green and red of "Landsat 8", when it is displayed on Erdas Imagine 2015 it appears as if there is a semi-transparent grey layer over it, the raster itself has no problems as I displayed it on ArcMap and it appears normal.