texts: Write RasterLayer by loop preserving the layer names and grouping models

texts:
I have a RasterStack object from witch I need to take some layers into a new RasterStack. To to that I've created a loop. My problem in that I am losing my original layers names in the process. As you will see below, when I take just one of the layers without the loop, the original name are kept.

In addition  I also need to group the layers (variables) into models, so I can iterate though the variables inside each model with a bigger loop I'm building.

How do I extract layers from a raster keeping the original layer naming?

Here is some exemple..

## stk1 is the original RasterStack
## stk1 has 3 models in it (.1,.2,.3) each model has 4 variables (day1, day2, day3, day 4)
names(stk1)
[1] "day1.1" "day2.1" "day3.1" "day4.1"
[5] "day1.2" "day2.2" "day3.2" "day4.2"
[9] "day1.3" "day2.3" "day3.3" "day4.3"

Here is the loop i've got running.

## Selecting day2 and day4 from stk1
variables <- as.factor(c("day2.1", "day4.1",
                         "day2.2", "day4.2",
                         "day2.3", "day4.3"))

for (i in 1:length(variables))
{
  writeRaster(stk1[[i]], filename = paste0("./data/stk1-", variables[i], ".grd"), format = "raster", bylayer = TRUE) # I've tried with and without bylayer.
}

stk1_select  <- stack(list.files("./data/", pattern = ".grd$", full.names = TRUE))
names(stk1_select) **# this names should be exactly as it is in variables**
[1] "day3.1" "day1.1"
[5] "day2.1" "day4.1"
[9] "day4.2" "day2.2"

Extracting one by one does keep the original layer name.

writeRaster(stk1[["day2.1"]], filename = paste0("./data/stk-", "day2.1", ".grd"), format = "raster")
day2.1 <- raster("./data/stk-day2.1.grd")
names(day2.1)
[1] "day2.1"

Finally what I need is like an array of rasters. stk1_select should be modelA( var: day2, day4), modelB (var: da2, day4), modelC (day2, day4)


texts: Dissolving shapefile python script not working, syntax help please

texts:
I do not understand why this simple script won't work. I can get it to dissolve in model builder but not from a script that I run from inside the map document. My output results are 4 polygons, just like my input file. I only want one.

Any suggestions?

# Import arcpy module
import arcpy
from arcpy import env

# Local variables:
env.workspace="C:\\Test\\shapefiles"

infeatures = "Multipolys.shp"
outfeatures = "C:\\Test\\shapefiles\\Dissolve.shp"

# Process: Dissolve
arcpy.Dissolve_management(infeatures, outfeatures, "", "", "SINGLE_PART", "DISSOLVE_LINES")

correct syntax is:

arcpy.Dissolve_management(infeatures, outfeatures, "", "", "MULTI_PART", "DISSOLVE_LINES")

ArcGIS 10.1 advanced.


texts: Updating MXD data paths on new computer: a few questions

texts:
Fairly new to python, trying to update the data paths for several dozen MXD files to a new computer that doesn't have the "B" drive of the author.

Does the physical\logical data drive have to actually exist on the new computer, for python to reassign using "mxd.findAndReplaceWorkspacePaths" (say if the MXD references a B: drive, do I need to worry that my data is on the F: or C: drive?

I have a fairly simple script, that I keep getting syntax errors on:

import arcpy
mapdoc = arcpy.mapping.MapDocument("CURRENT")

mapdoc.findAndReplaceWorkspacePaths("<the directory you are changing from>" , "<the directory you are changing to>")

mapdoc.save()

del mapdoc

Using: ArcGIS 10.3.1

Pseudo paths look like:

 from directory: "B:\Data_Directory\My Maps & Data\GIS_GDB_Layers\"

 to directory "C:\GIS_GDB_Layers\"

error: Parsing error SyntaxError: invalid syntax (line 1)

I don't have a B:\ drive on my computer- apparently the person who forwarded the data & MXDs did.


texts: Appending Feature Class from Geodatbase to another using ArcPy?

texts:
I'm trying to write a python script to load feature class from an existing geodatabase to a new one, here is my script:

import arcpy
from arcpy import env
import os
import traceback

#Set up workspace
outworkspace = r"C:\user\Geodatabases\Gdb_mena_v3.gdb"
workspace = r"C:\user\Geodatabases\gdb_4.3.0_v2.gdb"

arcpy.env.workspace = workspace

#get a list of feature classes
featureclasses = arcpy.ListFeatureClasses()

for fc in featureclasses:

    shp_file = fc

    if shp_file == "MAP":
    arcpy.Append_management(os.path.join(workspace, shp_file), "C:\\user\\Geodatabases\\Gdb_mena_v3.gdb\\MAP", "NO_TEST", "","")

    print('Your load is done')

    else: print('Error reading data')
print "Finished!"

when I run the script I get this message:

 Error reading data
 Error reading data
 Error reading data
 Error reading data
 Error reading data
 Your load is done
 Error reading data
 Error shape file did not get read
 Finished!

Note that in my first geodatabase I have 7 features classes


texts: Adding a PostgreSQL database to a web map

texts: I am interested in learning how to add a PostgreSQL database to a web map so I can do functions like select by attributes or locations. I am using Leaflet API to make the web map. What are the necessary steps to integrate a PostgreSQL database into a Leaflet?


texts: Mosaicking a selection from a raster catalog and returning the mosaicked raster

texts:
I have pretty much given up on trying to get help from ESRI forums so I thought I would try my luck with the good people here...

Here is what I have (in C#)...  I have an IFeatureCursor that is giving me a selection from an IRasterCatalog (so basically it is a cursor full of IRasterCatalogItems).  What I would like to do is take these rasters and mosaic them inside a method and return the mosaicked raster.

I have tried what I thought would be obvious beginning with an IRasterCollection...

IMosaicRaster m_mosaic = (IMosaicRaster)m_rasColl;
m_mosaic.WhereClause = "";

ISaveAs m_saveAS = (ISaveAs)m_mosaic;
IRasterDataset rasDataset =
       (IRasterDataset)m_saveAS.SaveAs("tmp_" + selectedFeatureClass.FeatureClassID.ToString(),
                                       tmpWS, "TIFF");

IRaster extractedRaster = rasDataset.CreateDefaultRaster();

// Compute the statistics for the raster.
IRasterCalcStatsHistogram m_calcStats = new RasterCalcStatsHistogramClass();
IStatsHistogram m_stHisto = new StatsHistogramClass();

m_calcStats.ComputeFromRaster(extractedRaster, 1, m_stHisto);

return extractedRaster;

This gave me what seemed to be a properly mosaicked raster, and the values seemed to be correct when I used the identify tool in ArcMap, however, the display looked like a raster that has not had statistics calculated and I can't seem to figure out why.  So I started looking at less obvious ways to do it but I haven't had much luck so far.

Suggestions?  Thanks a bunch!


texts: PostGIS - ST_PointOnSurface and ST_Centroid dropping Z

texts:
I am using Postgres 9.5.2 and PostGIS 2.2.2.

According to the documentation ST_PointOnSurface should not drop Z. This is what happens. If I run

SELECT ST_AsText(cell) AS pt FROM cells LIMIT 1

I get

POLYGON Z ((157.215247156568 3779.5275740661 76.0620819960913,157.215247156568 3937.00787401575 76.0620819960913,3.4336355447806e-07 3937.00787401575 66.9291338782187,3.43363565136201e-07 3779.5275740661 66.9291338782187,157.215247156568 3779.5275740661 76.0620819960913))

as expected (Z conserved). However, if I run

SELECT ST_AsText(ST_PointOnSurface(cell)) AS pt FROM cells LIMIT 1

or

SELECT ST_AsText(ST_Centroid(cell)) AS pt FROM cells LIMIT 1

I get

POINT(78.6076237499659 3858.26772404092)

and the 3rd coordinate is gone.

Am I doing something wrong?


texts: Buffer line between two coordinates with meters

texts:
I am trying to build a polygon (with long/lat coordinates) over a line between two points (long/lat as well). Basically what I used is line buffering with meters as an distance.
However I have noticed that my polygon is not exactly created with choosen buffer offset (e.g 100km as buffer gives me like ~75km range polygon on OpenStreetMap).
How can I achieve correct precision of buffer in meters?
I am using Java GeoTools.
Line builder:
public GeoLineBuilder() throws SchemaException {
    this.lineSchema = DataUtilities.createType(
            "",
            "Location",
            "locations:LineString:srid=4326," + "id:Integer"
    );
}

SimpleFeature buildLine(List<LocationPoint> points) {
    SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(lineSchema);
    GeometryFactory geometryFactory = new GeometryFactory();
    List<Coordinate> coordinates = points.stream().map(point -> new Coordinate(point.getLongtitude(), point.getLatitude())).collect(Collectors.toList());
    LineString line = geometryFactory.createLineString(coordinates.toArray(new Coordinate[] {}));
    featureBuilder.add(line);
    return featureBuilder.buildFeature(null);
}

Polygon over the line with buffer in meters:
private SimpleFeature bufferLine(SimpleFeature feature, long distance) throws Exception {
    GeometryAttribute gProp = feature.getDefaultGeometryProperty();
    CoordinateReferenceSystem origCRS = gProp.getDescriptor()
            .getCoordinateReferenceSystem();

    Geometry geom = (Geometry) feature.getDefaultGeometry();
    Geometry pGeom = geom;
    MathTransform toTransform, fromTransform = null;

    if (!(origCRS instanceof ProjectedCRS)) {

        Point c = geom.getCentroid();
        double x = c.getCoordinate().x;
        double y = c.getCoordinate().y;

        String code = "EPSG:3857";
        CoordinateReferenceSystem auto = CRS.decode(code, true);
        toTransform = CRS.findMathTransform(
                DefaultGeographicCRS.WGS84, auto);
        fromTransform = CRS.findMathTransform(auto,
                DefaultGeographicCRS.WGS84);
        pGeom = JTS.transform(geom, toTransform);
    }

    Geometry out = pGeom.buffer(distance);
    Geometry retGeom = out;

    if (!(origCRS instanceof ProjectedCRS)) {
            retGeom = JTS.transform(out, fromTransform);
    }

    SimpleFeatureType schema = feature.getFeatureType();
    SimpleFeatureTypeBuilder ftBuilder = new SimpleFeatureTypeBuilder();
    ftBuilder.setCRS(origCRS);
    ftBuilder.addAll(schema.getAttributeDescriptors());
    ftBuilder.setName(schema.getName());

    SimpleFeatureType nSchema = ftBuilder.buildFeatureType();
    SimpleFeatureBuilder builder = new SimpleFeatureBuilder(nSchema);
    List<Object> atts = feature.getAttributes();
    for (int i = 0; i < atts.size(); i++) {
        if (atts.get(i) instanceof Geometry) {
            atts.set(i, retGeom);
        }
    }
    return builder.buildFeature(null, atts.toArray());
}

Basically:

I am building a line in EPSG:4326 to have a line with long/lat coordinates
I am doing an EPSG:4326 -> EPSG:3857 transformation to be able to buffer by meters
Actual buffering
Reverse transformation to coordinates EPSG:4326 <- EPSG:3857
Displaying received polygon on the openstreetmap map

I saw that when I choose Poznan and Paris coordinates as a line pair - received polygon was for sure smaller than selected buffering (100000 meters), it was like ~75000 meters instead.
How to approach this problem?


texts: Using WK ID's when using CreateSpatialReference tool

texts:
I recently found out how to create specific spatial reference objects and add them to my arcgis python standalone script. sp=arcpy.CreateSpatialReference_management("GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],VERTCS['NAVD_1988',VDATUM['North_American_Vertical_Datum_1988'],PARAMETER['Vertical_Shift',0.0],PARAMETER['Direction',1.0],UNIT['Meter',1.0]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision","#","#","#","#","#","0")

As you can see, this is painfully long and is subject to errors. I would like to know how to go about shortening the input required to make this tool functional. Where can I place well-known ESRI ID's to simplify the inputs.


texts: Finding if a vertex is located between two edges

texts: I have a data set of polygons which I want to convert to lines in order to measure the length of each line. However, some of the polygons have vertices which are located not in a corner of the shape. Is there a way to determine if a vertex sits in a corner or not? I've tried Dissolve or Simplify without success.


texts: MapInfo "Error during Read: XML Library Load Error Required XML version unavailable"

texts:
I am using MapExtreme 7.0, I am running into the following error on an asp.net website (.net 3.5) running on a server with windows 2003 32bit.

"Error during Read: XML Library Load Error Required XML version unavailable" when trying to access some feature data.

I managed to track down some obscure reference in the 6.7.1 release notes and the fix was to install the latest "XML", I'm assuming that refers to MSXML, I have checked and the server has MSXML 6 installed.

Is there something else this may be referring to? Or is there another cause of the problem.

I managed to reflect the MapInfo dlls and have found the resource that contains the string but I cannot find where the error actually originates from


texts: I have a stack of MODIS images consisting of 400+ bands. A stack of NaN masks is created corresponding to each layer of the MODIS stack. I wonder if there is a way to apply the mask layers in bulk in ENVI on the MODIS stack instead of selecting just one mask layer to apply each time manually? If ENVI is limited, is there a way to bulk apply in other softwares/programs? I have access to ArcGIS, Matlab and R.


texts: Merging overlapping lines into one line. QGIS

texts:
I have a shapefile with multiple perfectly overlapping lines. Each line represents traffic at a certain point in time. I would like my output to have only one line but with the attributes from all other lines. Each segment has 123 overlapping lines

Any hints on how I should approach this?


texts: Examples of QgsFieldComboBox.setFilters() in building QGIS Plugin

texts:
I am making a plugin in QGIS which involves the use of QgsFieldComboBox, what I want to do is that the combo box only shows fields with data type of string (or to show fields with data type of integer, etc). Does anyone have any code example?

comboBox.setLayer(layer)
comboBox.setFilters(???)

Also, setFilters() takes QgsFieldProxyModel filters as parameter but when I do:

from qgis.core import QgsFieldProxyModel

It gives me an error saying that ImportError: cannot import name QgsFieldProxyModel


texts: Is it possible to use all Google Maps feature with Leaflet's GoogleMutant API

texts:
I'm using leaflet in my project for longer period of time.  Now it is required to use Street View feature of google map with it's Peg Men.

Is it feasible?

I've another question in mind where will I be allowed to consume all the feature google map's provide via Leaflet GoogleMutant API?


texts: updating field based on another field with arcpy.da.updateCursor

texts:
My code is very simple and it should work, but it doesn't. Perhaps another set of eyes can solve my issue...
I created new field (this part works), then look up value of another field and depending on it, update field "Surface" with either 3 or 2 (in my case value stays "null", but no error).

here the code:

import arcpy
from arcpy import env
arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = True

field5 = "SURFACE"

arcpy.AddField_management(input_fc, field5, "SHORT", "", "", "", "", "","", "Dm_VehicleCapacity")

fields2 = ['VehicleCarryingCapacityCurrent', 'SURFACE']
with arcpy.da.UpdateCursor(input_fc, fields2) as cursor:
    for row in cursor:
        if fields2[0] == 9:
            fields2[1] = 3
            cursor.updateRow(row)
        else:
            fields2[1] = 2
            cursor.updateRow(row)


texts: How to get a message when the map is moved in QGIS?

texts:
I want to get a message when every time I move the map.

I had:

def info(self):
    QMessageBox.information( self, self.tr( "Info" ), self.tr("info") )

I suppose that used to be something like that...

iface.mapCanvas().???.connect(self.info)


texts: arcpy join csv files to raster attribute table

texts:
In arcpy, I'm trying to join an list of .csv files with data of emissions per country, with the attribute table of a raster in order to make lookup tables later on. So far this is my script:

ras_countries = "E:\\eco_services\\data\\EU_28"

arcpy.env.workspace = "G:\\AIR_QUALITY\\paper\\GAINS_emission\\emission_gains"
landuse_tables = arcpy.ListTables ("emiss*")
for landuse in landuse_tables:
    tableViewName = landuse [6:-4]
    arcpy.MakeTableView_management(landuse, tableViewName)
    arcpy.AddJoin_management(ras_countries, "cou_abb", tableViewName, "cou_abb")

But I always get the same error message:

ERROR 000840: The value is not a Table View.
ERROR 000825: The value is not a layer or table view
ERROR 000840: The value is not a Raster Catalog Layer.
ERROR 000840: The value is not a Mosaic Layer.
WARNING 000970: The join field COU_ABB in the join table EU_28 is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.

I also tried converting the .csv files to .dbf but still no success.

any suggestions?

Failed to execute (AddJoin).


texts: Cartodb - Error: sql is read-only in NamedMaps

texts: My map is not a named map. Map privacy is set to only people with the link. Why I I getting this error when trying to do queries? Has anyone else run into this?


texts: GDAL error in QGIS: "Interpreting GDAL error: "gdal_rasterize: symbol lookup error: gdal_rasterize: undefined symbol: CheckExtensionConsistency"?

texts:
I'm using QGIS 3.4.6 with GDAL/OGR 2.3.1 to rasterize a points shapefile with the gdal_rasterize tool, and I keep getting the error 'Interpreting GDAL error: "gdal_rasterize: symbol lookup error: gdal_rasterize: undefined symbol: CheckExtensionConsistency"'.

What's going on here?  How do I troubleshoot this?


texts: How to determine the status of snapping in QGIS 3?

texts: I can use toggleEnabled to turn snapping on and off, but is there anyway to know the status of the snapping tool (on or off)? as toggleEnabled doesn't return a value nor does it take a value.


texts:
I am trying to add a CartoDB basemap to QGIS and am 1) trying to figure out if it is possible and 2) trying to figure out how to do it.

I've installed the QGIS CartoDB plug-in and have connected to my account. When I select Add CartoDB layer from Web >> CartoDB Plugin >> Add CartoDB Layer I only see the data layers I've uploaded to CartoDB, not any base maps. How can I access the basemaps?

Even after installing QuickMapServices I do not see that CartoDB base map:


texts: Labels on the QGIS WMS layers are not drawn according to the style definition

texts:
I have an issue regarding the styling of labels of WMS layers using QGIS server. I can't seem to get them to display how I want.

What I did:

create a QGIS project (2.2 on Mac) with a layer from a database
style the layer, including the labels on QGIS desktop
save project on the server (Ubuntu) to create the WMS

What I get:

my WMS layer, nicely showing up
except for the labels that seem to use a default symbolization, whatever the changes I make in the project (understand big black labels)

What I wished:

control the size and maybe also the font of the labels.

Is there a nifty parameter that I forgot somewhere or is the labeling engine not working propoperly with WMS? Any hints or similar experiences welcome.


texts: Cannot Edit Python Script from ArcMap/ArcCatalog with IDLE?

texts:
I'm trying to set my Script Tool Editor in the Geoprocessing Options dialog so that I can edit a script tool within ArcMap. Currently, whenever I right-click on my script and select "Edit..." from the resulting menu, nothing happens. I've set my editor to what I think is the correct .exe file (C:\Python27\ArcGIS10.1\Lib\idlelib\idle.pyw), but selecting "Edit..." on the script tool menu still has no effect.

My issue seems to be similar to Cannot Edit Python Script tool using its context menu from ArcMap/ArcCatalog? which was never answered but I am using ArcMap 10.1 SP for Desktop Basic and Python 2.7 (I think).


texts: How to prevent Leaflet from loading tiles from server

texts:
I am currently developing a couple of solutions with LeafletJS and it frustrates me that my Firefox browser loads the tiles from the server (e.g. Stamen or OSM) and not from cache every time I reload the page (I only reload the page to run my latest js code, not to see fresh tiles).

I used to work with OpenLayers and didn't have that issue - on the contrary, I added some extra code to force OL to get fresh tiles in specific cases.

I have Leaflet version 0.7.3 and Firefox version 38.0 on Linux Mint Mate 17.1 64 bit.


texts:
I have two polyline shapefiles, Green and Gold which are supposed to match but don't. Gold is larger than green, and slightly offset. Some line IDs are a match, but other IDs have been changed.

I am trying to combine these in ArcGIS with the following steps:

"Join" Gold to Green (Green target) on ID = segment_id
Select all Green segments where segment_id is not null and create "matched_segments" layer
Select all Green segments where segment_id is null and create "unmatched_segments" layer; delete joined attributes
"Spatial Join" Gold to "unmatched_segments" layer to create "unmatched_spatial_map" layer
Combine "unmatched_spatial_map" and "matched_segments" into one layer and export

At step 4 (completed) I have two layers where each segment is contained in either one layer or the other (not shared) and would like to combine them to one layer with Append (from Data Management Tools). This is not working. "pipes_matched" has the prefix 'gold' for all the field names which were joined from the input file, and the field name is greyed out and unable to edit. When I run Append with 'Input fields must match target fields' I get ERROR 000466 for schemas that don't match. When I run with 'Use the field map to reconcile field differences' I also get an error.
How can I solve this so that I have a single layer (to export to a shape file) which is joined with the gold GIS data on ID when there is a map and on the closest pipe when there isn't?


texts: How do I write a MYSQL spatial query to find all polygons (rectangle) that contain any points in a table?

texts:
Something like:

SELECT * from polygons WHERE bounds CONTAINS(SELECT geom_point from points)

Can someone help me out with the syntax? Bounds will always be a rectangular polygon and everything is in lon/lat. I could just check that each points lat/lng is between the min/max lat/lng of each polygon - but I cannot figure out the syntax to get MYSQL to loop / check all polys vs all points.

Thanks in advance.


texts:
I am new to GIS and I need some help with the architecture of a new system.
The UI will be written using javascript with openlayers 3.

The server should be written in JAVA and as I understand, it can't hold the maps on it's own (too complicated), but has to query and update some kind of a maps server. I'd like some help choosing the right maps server.

A friend recommended mapguide - but I see it's mapguide4j is too complicated and it uses Play framework - We want to use Spring.
Or maybe I didn't get that correctly?


texts: Display data from postgresql on click map

texts:
I work with openlayers and geoserver. I'm searching for code that when I click on a point of a layer on the map, a popup window appears that contains information from the database (postgres) for this point.

Can you help me please?


texts: Bezier curves imported with 64bit MAPublisher become angular?

texts:
The bezier curves in ArcGIS (10.4.1) become angular when I import them with 64bit MAPublisher 9.8 in Illustrator (CC2015)

Are there any settings to fix this?

.


texts:
I have huge polygons (expanded on more than 25 degrees) with vertices positioned on degrees coordinates grid nodes on WGS84 CRS :

I'd like to calculate, in the most exact way, the area of this polygon. In my case I use the UTM42S CRS.

In QGIS, I convert the polygon from WGS84 to UTM42S CRS
Before I calculate the area, I display the degrees coordinates grid layer (WGS84) and my new polygon layer (UTM42S) with a project CRS set on UTM42S with the "on fly" option. We can see that the longer polygon line doesn't follow the grid line. It means that the area will be greater than the reality.

The only solution I have is to add other polygon vertices on grid nodes (when it's possible) on WGS84. Like that after the conversion to UTM42S strait lines will more follow the grid curve lines. Any ideas ?

Here are the polygons area into UM42S CRS :

FID area-1
0   9 213
1   2 692
2   8 271
3   39 649
4   273 653
5   55 684


texts: Using a shapefile as mask with rasterio mask gives different affine transforms with rasters of the same projection and resolution

texts:
I have a list of n file paths to geotiffs of the same resolution and projection but with varying extents. In brief the goal is to

use Rasterio's mask module to use a shapefile containing a single polygon feature that is the mutual overlap of all the geotiffs to open each geotiff cropped to the the mask as an array,
use numpy vstack to build a 3d array with a depth of n, and
use the the affine transform of the masked array to write out a multiband raster.

The issue I'm running into is that mask function returns out different affine transforms for thedifferent geotiffs. Given the same proj and resolution and same mask geometry, I expected that the affine transform returned from the mask function would be the same and that could be used to save out the multiband raster. Since that is not the case, using the a single affine transform to write the arrays back to raster results in misalignment between the bands. When just looking at the source geotiffs in QGIS they all seem aligned. Here's some code of the approach, example output of the different transforms, and some of the source raster info
import geopandas as gpd
from rasterio.mask import mask

def geofunc(listofpaths):
    arrlist = []
    afflist = []
    for pth in listofpaths:
        arr, aff = mask(rio.open(pth), shapes=msk.geometry, crop=True)
        arrlist.append(arr)
        afflist.append(aff)
    return arrlist, afflist

msk = gpd.read_file(aoishp)  # aoishp is a path to a shapefile
arrs, affs = geofunc([path1, path2, path_n])

## The mismatched transforms
affs[0]
Affine(12.5, 0.0, 328675.256933,
       0.0, -12.5, 6822881.355124)
affs[1]
Affine(12.5, 0.0, 328669.934103,
       0.0, -12.5, 6822872.680881)

## Some info of the rasters
rio.open(proclist[0]).crs       #  proclist is the list of paths
CRS.from_epsg(32612)
rio.open(proclist[1]).crs
CRS.from_epsg(32612)
rio.open(proclist[1]).bounds
BoundingBox(left=328569.934103, bottom=6776610.180881, right=419257.434103, top=6851810.180881)
rio.open(proclist[0]).bounds
BoundingBox(left=328675.256933, bottom=6776693.855124, right=419112.756933, top=6851706.355124)
rio.open(proclist[0]).res
(12.5, 12.5)
rio.open(proclist[1]).res
(12.5, 12.5)


texts: Determining the problem with an WMS layer in Geoserver

texts:
I would like to ask a question as I'm having an issue here which I cannot really determine what is the actual problem.

Is it that having more styles per layer is not really possible or could it be something else?

Oddly enough, I noticed that if I have more than 6 styles the others don't show up. However, when I looked at the Start Geoserver window there was an error that caught my attention: "The requested style cannot be used with this layer. This style specifies attribute of p_T20 and the layer is..." and it goes on with  multiple errors being mentioned.

For this exercise, I'm running on a local server on my machine.


texts: How to create a constant numeration of ID (indexing) in an attribute table (ArcMap)

texts:
Is it possible to create a constant indexing (as a numbers) in an attribute table of the theme with many objects that manner that it will remain (will not change) during update of this theme?
I found the instrument “Drop Index (Coverage)”, but it seems that it works only with ArcInfo license. Is it possible to create “constant” numeration of ID in ArcMap?


texts: How can I get GeoJSON from Esri software to load in Leaflet?

texts:
So I am trying to take a GeoJSON file that I have generated from an Esri shapefile in ArcGIS Pro, but it is not rendering on the map for some reason...

The GeoJSON that was created from ArcGIS Pro is here. I have included my current JavaScript below as well. I feel like there might be something wrong with the way the data is formatted or something, I don't really know, as I am unfamiliar with all of this. I left in some commented code that didn't seem to do anything differently, but maybe it can work if a small change is made.

$(document).ready(function () {
    let map = L.map("map").setView([19.042805, -70.581183], 7.2);

    L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
        attribution: '&copy; OpenStreetMap contributors'
    }).addTo(map);

    $.ajax({
        dataType: "json",
        url: "/static/yaque_del_norte_viewer/GeoJSON/yaque_del_norte.geojson",
        success: function (data) {
            console.log(data);

            let myStyle = {
                "color": "#2427ff",
                "weight": 5,
                "opacity": 0.65
            };

            L.geoJSON(data, {style: myStyle}).addTo(map);

            // let yaqueDelNorteLayer = new L.GeoJSON();
            //
            // $(data.features).each(function (key, data) {
            //
            //     console.log(data);
            //     yaqueDelNorteLayer.addData(data);
            // });

            // console.log(yaqueDelNorteLayer);
            //

            // yaqueDelNorteLayer.setStyle(myStyle);
            //
            // yaqueDelNorteLayer.addTo(map);
        }
    });
});


texts: How to perform euclidean distance analysis with coefficients?

texts:
I would like to perform an euclidean distance analysis in ArcGIS 10 for multi-points which represent seaports. The thing is these facilities have different capabilities for handling and i want to take this specialty into consideration in my analysis. Thats why i added a column in my attribute table.

How i represent this coefficient in my analysis?

Thanks for your helps.


texts: Will Esri’s Road Closures solution perform a reverse geocode to get Waze street names?

texts:
We are new to the Waze for Cities program. Currently, I am trying to get information regarding how to submit road closure data to Waze using Road Closure solution. I read all the documentation for this solution, as well as the guidelines provided by Waze for this task. Being an Esri shop, we want to use this tool rather than Waze Map Editor since WME will not work for our workflow.

I noticed that Waze Map Editor tool will auto-populate the street name information when creating road closures. Since we are using the Esri solution tool for Waze, that is, the web mapping application for road closures, street names are not auto populated. As you may know, the streets shown in the Esri’s base map may vary from Waze basemap street names. Will Esri’s solution perform a reverse geocode operation using as an input the geometry delineated in the web map to get Waze street names?

See here: https://developers.google.com/waze/data-feed/problem-resolution#reverse-geocoding-api


texts: How do I add all features of a Geometric Network to the the dataframe with one parameter in a custom tool that uses python?

texts: I want to add all the Feature Classes of a Geometric Network to my data frame with python.  I want the user to be able to select the dataset of their choosing in a tool, so it would use one parameter.  I want to then use the added Feature Classes to the data frame to run the rest of the script, which I am pretty sure I can figure out once I get the layers in the map.  Does anyone have any ideas? Sorry, this is probably a rudimentary question, but I am new to python.


texts: Does ArcGIS Online allow Select By Location?

texts:
This question refers to the standard ArcGIS Online map, eg http://your-site.maps.arcgis.com/home/webmap/viewer.html

Is it possible to perform a Select By Location, eg to find points from one layer which fall within a polygon layer?

There are a bunch of operations available under the Analysis menu but I can't see Select By Location. Is it necessary to build something custom (eg using the ArcGIS JS API, Web AppBuilder, etc) in order to get the Select By Location functionality?


texts: ArcPy block statistics need to output different resolution

texts:
I am attempting to script a BlockStatistics where the outRaster has a resolution different from the input raster. Input raster cells are approx 5.8 cm, and I want to resample taking the mean of a 30x30 window, returning approx resolution 1.7m x 1.7m. I can do this in ArcGIS Pro with BlockStatistics tool and changing environment to cell size to 1.772 but for the life of me cannot find where to change the environment for the tool using ArcPy. Also, if there is a simpler solution I am missing, I am open to it.
OutRas returned with the following code snippet calculates block statistics fine, but returns raster of 30x30 blocks with the identical mean calculations and original resolution.
OutRas = arcpy.sa.BlockStatistics(inRas, "Rectangle 30 30 Cell", "MEAN", "")


texts: Create cross section with polygon intersects in QGIS

texts:
Is it possible to create a cross section that shows the intersected polygons on the relief (in QGIS 2.0.1)? I would like to try this out for a geologic map.

I'm using VoGisProfilTool:
http://plugins.qgis.org/plugins/VoGisProfilTool/
and a NASA-DEM


texts: How to group different values of a variable so that they plot with same shape?

texts:
I have a table with two columns of number X and Y, and one column of strings called crystals.

crystals
 A
 A
 B
 C
 B
 D
 B
 C

Using ggplot2, how can I make a scatter plot of X versus Y with two types of points only, one shape for Crystals A, B, C and one shape for crystals D?

Tried to create subset but have not found a way to make 2 groups out of the 4 different values.
I also tried to assign manually the same shape for A, B and C but then get a legend with the 4 entries instead of 2.


texts: Styling a WFS layer, how to control title in legend?

texts:
I have a problem with styling a WFS in OpenLayers. I tried different ways (from the examples) to style the WFS, but I could only get the default style in orange. Then I found one way, which worked for me. The problem is that in the legend is written "Untitled 1" next to the now correct symbol as with the default style. How can I change this title or remove it?

Has anybody have a solution on how to style the WFS with an external SLD? I also tried different examples found on the web but I couldn't get it to work.

Here's the code I am using to style the WFS but that displays the wrong title:

var StyleObjekt = new OpenLayers.StyleMap({
                            "default": new OpenLayers.Style({
                                                          fillColor: "#fff000",
                                                          strokeColor: "#666666",
                            }),
                            "select": new OpenLayers.Style({
                                                            fillColor: "#ffffff",
                                                            strokeColor: "#666666",
                                                            }),
                            "temporary": new OpenLayers.Style({
                                                            fillColor: "#ffffff",
                                                            strokeColor: "#666666",
                                                            })
            });

 var gwg_f_wfs = new OpenLayers.Layer.Vector("Flurstuecke_WFS", {
            styleMap: StyleObjekt,
            strategies: [new OpenLayers.Strategy.BBOX()],
            protocol: new OpenLayers.Protocol.WFS({
                version: "1.1.0",
                url: "http://localhost:8080/geoserver/wfs",
                featurePrefix: "NewGeo",
                featureType: "Flurstuecke_GWG",
                featureNS: "http://www.newgeo.de/GeoServer",
                srsName: "EPSG:2398"
            })
        });


texts: How to standardise colour-value graduations across multiple maps

texts:
In QGIS I have a vector dataset (imported CSV vector + csvt) spanning multiple years, and a shapefile shape which is essentially a grid.  By joining the CSV file with the shapefile I can apply a colour to each grid square and hence a gradient scheme for the whole grid. I have created a series of maps - one map per year, each with its own colour gradient.  However the colour gradients are not comparable between maps because the range of values is different among years. How can I make the gradient scales standardised across all maps?

I've tried splitting the CSV file into years, and with each year individually, it does not seem possible to create a 'custom' range of values for the gradient, say 0-100, even though the values span say 20-70. I've also tried working with one CSV file containing all years, thinking that I could somehow 'trick' QGIS into applying a gradient that spans the range of the whole dataset, while at the same time only map individual years.  Categories doesn't seem to be an option as if I have 200 grid squares I get 200 categories each with own unique value.

Any help would be much appreciated.  It may be that there is a database-based solution which I am open to.  I am also open to GRASS solutions, but I am less familiar with GRASS as this stage.


texts: Can someone explain to me how to use wild cards with QuickOSM? I'm trying to match all shops with a name starting with a specific word (e.g. Costco).


texts: QGIS 3.0 Field calculator - Call second layer in Intersect

texts:
I am trying to use the 'intersects' function within an 'if' condition. If a polygon in Layer A intersects a polygon in layer B, then the cell should display 'YES', otherwise 'No'.

However, I don't know how to properly call the Layer B in the field calculator. The expression I am working on looks like this :

if(intersects( $geometry,geometry('LayerB') ),"YES","NO")

Does anyone know how to put it correctly?


texts: Unable to expand projected co-ordinate system folder in ArcMap 10.1

texts:
On trying to re-project any raster from GCS to any PCS, in the "Spatial Reference Properties" window, my "Projected Co-ordinate Systems" folder is not expanding. Where as the GSC folder is expanding. Even I tried to clear the spatial property of the layer and tried to set it to the desired PCS, there also its not expanding.

After restarting ArcMap, before clicking on the ArcToolBox, if I go to layer properties, there it expands. Once I clicled on ArcToolBox, it doesn;t expand either from Layer properties or from Project Raster Tool. I have attached the the image of my window. If at all any one can guide me a solution. Please.. please.. please... I'm stuck almost since last 3 months.


texts: BingService doesn't accept requests for Hybrid and Other Bing Map formats

texts:
We have used BingService with the help of following code:

String baseURL = "http://ak.dynamic.t1.tiles.virtualearth.net/comp/ch/${code}?mkt=de-de&it=G,VE,BX,L,LA&shading=hill&og=78&n=z"`;

map.layers().add(0, (new TileLayer(new BingService("Hybrid", baseURL))));

but BingService always create a Road map.

Is there any way to create Hybrid and Arial map?


texts: Coloring heat map, with custom characteristics

texts:
I need to color a "heat map", with these characteristics: "Pseudocolor monoband", "Discrete" and with two levels:

<= "STATISTICS_MEAN" -> Color 255/255/255
"Inf" -> Color 196/255/185

I have to do it in different "heat maps" and I want to do it with the "Graphic Modeler".
How can I do it?.


texts: How to resolve error in Orfeo Toolbox for QGIS 2.18 - "option -ram does not exist in the application."?

texts:
I followed the instructions to install the Orfeotoolbox for QGIS 2.18 on Windows, detailed in this post. Everything worked, and the appropriate algorithms are showing under the Processing Toolbox:

However, while some algorithms are running successfully, when running many of the tools, I get this error:

  ERROR: option -ram does not exist in the application.

So far, the error occurs in these algorithms:

Exact Large-Scale Mean-Shift segmentation, step 3
Compute Images second order statistics

I have tried re-installing QGIS and the toolbox on a different computer, and had the exact same problem. I've also tried clearing the parameter in the tool interface. Has anyone had this problem and/or knows a workaround?


texts: Which UTM system do you recommend for Tamaulipas, Mexico?

texts:
I'm working with QGIS 1.8.0. Could anyone tell me which UTM system would be better to project a map of Mexico? There are too many to choose from. I have to convert my maps from WGS84 to UTM to do some distance measurements.

Please let me know.
Regards!


texts: How to locate and list peaks of a specified height in OpenStreetMap?

texts:
Is it possible to calculate a list of peaks for a specific area given a specific prominence (say 200 meters)?

If so can this be achieved in an automatic, or semi-automatic, way?


texts:
I am working on a library project for rectifying old map imaging on a map. The images are in tiff and jpeg formats.

After rectifying the images I want to save the geographical information in a PostGIS database, while the image itself will remain in another DB, in the same format.

What methods I can use to query the images DB with geographic data based query? (i.e. polygon search on the map will retrieve all rectified maps on that area)


texts: Calculating street distances using ArcGIS Desktop?

texts:
I am trying to calculate the length of a street.

I have a census database (https://www.census.gov/geo/maps-data/data/tiger-line.html) of streets that I am viewing in ArcGIS, but the actual length of the streets is missing, although it must be somewhere if it is being displayed on the map.

Does anyone know of any API or tool that will take a street address (state and county included) and output the how long the road is?


texts: GeoExt printPage cannot handle WFS vector features

texts:
I recently added the print function to my GeoExt web application, using

<script type="text/javascript"src="http://localhost:8080/geoserver/pdf/info.json?var=printCapabilities"></script>

to define the printCapabilities variable,

        var printProvider = new GeoExt.data.PrintProvider({
            method: "GET",
            url: "localhost/geoserver/pdf",
            capabilities: printCapabilities
        });

for the print provider and within the mapPanel

bbar: ["->", {
                text: "Print...",
                handler: function() {
                    var printDialog = new Ext.Window({
                            items: [new GeoExt.PrintMapPanel({
                                sourceMap: mapPanel,
                                printProvider: printProvider,
                                    listeners: {
                                    "loadcapabilities": function() {
                                    var printPage = new GeoExt.data.PrintPage({
                                        printProvider: printProvider,
                                        customParams: {
                                            mapTitle: "LVR WebGIS",
                                            comment: "created with GeoExt"
                                        }
                                    });
                                        printPage.fit(mapPanel, true);
                                        printProvider.print(mapPanel, printPage);
                                    }
                                    },

                                })
                            ],
                         bbar: [{
                                text: "Create PDF",
                                handler: function() {
                                    printDialog.items.get(0).print();
                                }
                            }]
                    });
                printDialog.show();
                }
            }],

to add the function to my mapPanel. Regarding this, I have two issues to be solved:

1.) the printPage popping up does not contain anything but the button saying "Create PDF". There should be a preview of the map to print, right?

2.) when some WFS vector features are displayed in the map, the print function does not work properly anymore. It opens a new tab in my browser but does not create the PDF. Instead, the new tab is left about:blank.

What am I doing wrong here? Is the information I provided enough for you guys to help me solve the problem? Thanks a lot for your help!


texts: QGIS 2.14 Symbol Colors Changing on Zoom

texts:
I'm having this problem with Rule Based Symbols.
I am able to change the color of a dot based on the name field but when I zoom in and out the colors bounce around to other dots.
By zooming in and out of my map, the colors change with no regards to the color rule I set.
I can not figure it out, maybe its a rendering thing but I'm not sure. Attached are screenshots. All these dots are supposed to be the brown like color. If I zoom in closer they will all turn brown.


texts: Right table values missing from AddJoin() when Join table order reversed?

texts:
Trying to get to the bottom of an unusual behaviour which I am seeing in a join. Background as follows:

I have a feature class of centroids which have been located on the network.

I have a solved OD Cost Matrix layer.

I now want to join the OD lines sub-layer to the centroids so that I can access the centroid snap-distance.

The join works in one direction just fine and returns values from the left and right tables.

However when the order of the join is reveresed, while the join still returns results, only values from the left table are included but not the right. I am wondering as to why this is the case?
    # Create a temporary feature layer to allow us to add a join to the selected centroids.
    centroid_layer_temp = arcpy.MakeFeatureLayer_management(selected_centroids,'centroids_layer_temp')
    lyr_desc = arcpy.Describe(centroid_layer_temp)

    #(1) This join works and returns joined rows as expected. Values from both left and right tables are included.
    matched_centroids = arcpy.AddJoin_management(centroid_layer_temp, "OBJECTID", od_line_layer, "OriginID", "KEEP_COMMON")

    #(2) This join works i.e. it returns rows but only values from the left table are included in the result...why?
    matched_centroids = arcpy.AddJoin_management(od_line_layer, "OriginID", centroid_layer_temp, "OBJECTID", "KEEP_COMMON")

    # Get centroid data and store for later use.
    # This throws an error for (2) above because no fields from the right table are included in the result.
    grid_centroids = []
    centroid_search_fields = ['{}.OBJECTID'.format(lyr_desc.name), '{}.Distance'.format(lyr_desc.name)]
    with arcpy.da.SearchCursor(matched_centroids, centroid_search_fields) as search_cursor:
        for row in search_cursor:
            grid_centroids.append(row)

    # notify is a convenience method to log and output to UI at same time.
    Helpers.notify(grid_centroids)


texts: Using CRS EPSG:3857 but misalignment between Stamen background and coordinates of polygons in GeoPandas

texts:
Using GeoPandas and Contextily to create maps with a Stamen background map I notice a misalignment between the background map and the coordinates of polygons drawn on the map. The shapefile is from the official Dutch Bureau of Statistics so the coordinates are correct.

What can I do?

My usage should be inline with the instruction and the example notebook. Others seem to have the same problem, e.g. this tutorial.

My original projection is:

{'proj': 'sterea',
 'lat_0': 52.15616055555555,
 'lon_0': 5.38763888888889,
 'k': 0.9999079,
 'x_0': 155000,
 'y_0': 463000,
 'ellps': 'bessel',
 'units': 'm',
 'no_defs': True}

I perform:

dfmap = dfmap.to_crs(epsg=3857)

The crs becomes:

{'init': 'epsg:3857', 'no_defs': True}

I then perform:

ax = df.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
ctx.add_basemap(ax)

But the basemap and polygon are misaligned, as per the map image (the polygon should be more to the left and a bit lower so that the curve aligns with the water on the map):

I tried different zoomlevels, checking axis values, setting ax.axis('square'),...

On axis values, adding the basemap changes the values:

print(ax.axis())
ctx.add_basemap(ax, zoom=15)
print(ax.axis())

Returns:

(542505.2085419805, 543235.6744236776, 6866130.598686902, 6866937.041798945)
(541785.6564853293, 544231.6413904548, 6865879.628687669, 6867102.621140235)

Which I would not expect considering the following in add_basemap(ax)?:

xmin, xmax, ymin, ymax = ax.axis()
# restore original x/y limits
ax.axis((xmin, xmax, ymin, ymax))

I also tried to add the basemap after performing dfmap.to_crs(epsg=4326) but then I get the sea as background. And I also tried with ll=True as well as ll=False in the ctx.bounds2img function:

def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url, ll=True)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

I lack the knowledge about long lat vs xy but if I read about how "the Web Mercator uses the spherical formulas at all scales whereas large-scale Mercator maps normally use the ellipsoidal form of the projection. The discrepancy is imperceptible at the global scale but causes maps of local areas to deviate slightly from true ellipsoidal Mercator maps at the same scale. This deviation becomes more pronounced further from the equator, and can reach as much as 40 km on the ground." then I guess it relates to that.

But how?


texts: How to disable map navigation in the ArcGIS Server JavaScript API v4.x

texts:
In the older version 3.x of the ArcGIS Server JavaScript API, it was possible to disable map navigation using disableMapNavigation() as explained in Disabling map in ArcGIS API for Javascript
What is the equivalent for disabling map navigation in version 4.x of the ArcGIS Server JavaScript API? The use-case is to prevent map navigation temporarily, while an external query is running.
The mapView has a navigation property which allows certain types of navigation to be disabled, as in the example at https://developers.arcgis.com/javascript/latest/sample-code/sandbox/index.html?sample=view-disable-navigation, but this doesn't seem to allow all navigation to be disabled.
Eg, even when using:
navigation: {
  gamepad: {enabled: false},
  browserTouchPanEnabled: false,
  momentumEnabled: false,
  mouseWheelZoomEnabled: false
}

the user is able to pan the map with the mouse, zoom by double-clicking, rotate the map using the right-mouse-button, etc.


texts: Is it possible to display related tables information (1-M relationships) on the layout while using data driven pages?

texts:
I'm using data driven pages to generate a report for every point that I have (each point represents a store) and I have used Dynamic Text to display the tabular information of these points like store name and size... but I have some information in a table related to this feature by 1-M relationships like names of employees in each store and their ages.

So how can I display the related table information on the layout and make it dynamic so that when I change the page using the data-driven page the info will change also. does Production Mapping extension helps in this case and how?


texts: Add bands or new raster column in Postgis?

texts:
Would it be faster to work with a new band or a new raster column for a raster dataset in Postgres with new pixel values, assuming that the size, extent and everything else relevant remains constant?
For example: If I've reclassified the raster column and want to retain the original values - should it go into a new band or just another 'reclass_rast'column?
Concerned with querying capabilities and size on disk.


texts: MODIS resolution missmatch between MRT and GDAL processing

texts:
I've been using MRT (Modis Reprojectin tool) to extract EVI bands into a geotiff with geographic coordinates (epsg4326). I did not specify an output resolution to MRT. The output from MRT had a spatial resolution of 0.002543280392921 degrees.

Recently I tried using GDAL to extract EVI bands. First I used gdal_translate to extract the EVI band in SIN. That gets me a image with 231.65m resolution (native MODIS resolution - OK). But when I use gdalwarp to convert the image to epsg:4326, the output resolution is 0.001682583113332 degrees! (again, I did not specify the output resolution in gdalwarp). If I convert to UTM, it gives me a resolution of 179.7m.

Shouldn't gdalwarp give me the same resolution obtained from MRT?


texts: Loading single feature from server - layer property is missing

texts:
I have a layer that consists of a few 100 features, of which only those are loaded that fit into a bounding box. Later, I want to check on the server side, whether a specific feature with a specific attribute exists and then zoom to that feature.

I am not sure about the appropriate way how to do this, so I just call the read function of my layer's protocol, with a filter parameter. The result looks fine at first, but the retrieved feature seems to be inclompete, it lacks for example a vale for the layer property.

Do I need some other attributes on the read-functions' option parameter? Or is there a simpler approach to the problem?

var mylayer = new OpenLayers.Layer.Vector("WFS", {
        strategies : [new OpenLayers.Strategy.BBOX(), saveStrategy],
        protocol : new OpenLayers.Protocol.WFS({
            url : baseurl + "/wfs",
            version : "1.1.0",
            featureType : "myfeature",
            featureNS : 'myNS',
            srsName :"EPSG:4326",

        }),
        renderers : ["Canvas"],
    });

 ......

  var filter = new OpenLayers.Filter.Comparison({
    type : OpenLayers.Filter.Comparison.EQUAL_TO,
    property : "myProp",
    value : 123
    });

    mylayer.protocol.read({

      filter : filter,

      callback : function(resp) {

         console.log(resp.features[0]);

      }

    });


texts: Create KML "animation" / time domain from separate KML files

texts:
I'm currently working on creating a weather map system with leaflet. I'm using grads to create KML files from NOAA/NCEP GFS data in grib2 format. Grads seems to be only able to write out one layer pro KML, so I end up with 60 KML files - one for each time index.

Is there any Linux command line tool which could combine those 60 KML files into one and also write the correct metadata so I could play a time domain animation of the file in Leaflet or Google Earth?


texts: Extracting point data from Digital Raster Graphic (DRG) in ArcMap?

texts:
I am trying to run an IDW analysis based off stream and springs data.

I have the hydrography.eoo file from GeoComm but am now wanting to create a point layer of the springs in the area which show up on the DRG.

Is there any way to extract point data from a Digital Raster Graphic (DRG) in ArcMap without manually creating a spring layer?

I'm using ArcGIS 10.1 for Desktop.


texts: R conditional logic embedded within calc function for raster

texts:
I am trying to create a new raster from one existing raster using conditional logic, ifelse statements. The existing raster has numeric values and NA values. Here is what I tried so far, based on other posts within Stack Exchange:

f <- function(phos_all) {
        ifelse(is.na(phos_all) | NA,
          ifelse(phos_all<=4,30,
            ifelse(phos_all>=4) & (phos_all<10), 20,
              ifelse((phos_all>=10),0)))
}

phos_fert_2.5 <- calc(phos_all, fun = f, na.rm = TRUE)

When I run this code I get the following error:

  Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
    cannot use this function. Perhaps add '...' or 'na.rm' to the function arguments?


texts: What does it mean if a topographic map is indicated to be "West of fifth Meridian"

texts:
The Natural Resources Canada topographic map A 741, 83 C/7, edition 3 has a (bilingual) indication at the bottom stating:

  WEST OF FIFTH MERIDIAN       OUEST DU CINQIÈME MÉRIDIEN

What does this mean?  It's around 117°W, so it is indeed west of the 5th meridian, but I think something else is meant here.  Does this have a specific meaning on a topographic map?


texts: How do I request an output from a WMS service in a non-standard projection?

texts:
I want to use a map service which is available in WMS format. The server is Mapserver, and the spatialreference of the map service is 4326.

I wish to use it in a slippy map (web based Javascript Map), with a custom projection. This projection does not have a EPSG code. Is it possible to request the output, by providing the well known text for the projection?

Update

I can see that the GetMap request has a SRS paramter, and usually an EPSG code is supplied, like this: SRS=EPSG:4326. Can we pass something else to this SRS Paramter?


texts: Misalignment of EPSG:4326 and EPSG:3399 coordinates

texts:
I'm having a QGIS project in local coordinates (Gauss-Kruger, EPSG:3399) with some layers in the same reference system. On top of that, I want to display point coordinates, which I digitized using google maps, so I assume them to be in Pseudo-Mercator/ WGS84 (EPSG:3857).

Now, when I add the point coordinates, they won't show up in the map. If I change their reference to EPSG:4326, they are displayed, but they have a mismatch of more than 150 meters.

In addition, I loaded a WMS layer, in WGS84 (EPSG:4326) which fits fine with the point coordinates but also does not align with the other layers from the project.

I included a partial view of the project, where you can see the misalignment. The overlay shows a river obtained from a DTM (in local coordinates) in blue. The displacement with the river in the orthoimage is clearly visible. The green points are the ditigized ones, as said above.

I am aware, that the on-the-fly transformation in the GIS might not be able to get me the exact match between local coordinates and global WGS-coordinates. But somehow I need to get the local coordinates from the point locations without having to digitize them again, because I am also missing georeferenced orthoimages in the local coordinates.

Any ideas how to handle this?


texts: Leaflet: Scroll wheel controlling map zoom and not drop down menu

texts:
Edited 6/19/2015

I have a custom leaflet control box with a drop down menu.  When I use the scroll wheel over the menu, the map changes zoom level instead of scrolling the menu.  The problem occurs on IE and Firefox, but not Chrome.

How can I let the user operate the drop menus with their scroll wheel, when their cursor is over the menu?

        <link rel="stylesheet" href="http://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.3/leaflet.css" />


texts:
I have exported a webmap created in qgis2web (QGIS 2.18.13) open layers. There is a list of attributes that I wanted to include in the popup, which are mostly related to a country's demographic and socioeconomic profile. From the image below however, you will see that the popup is pretty big, out of range and without a scroll.

I have attempted to edit the qgis2web.css file within resources to try and address some of the points but I am completely new to CSS so am unfamiliar with the syntax... Could any of you please help to:

Adjust the width/height of the popup
Ability to scroll through the popup

html, body {
    height: 100%;
    padding: 0;
    margin: 0;
    font-family: sans-serif;
    font-size: small;
}

th, td {
    vertical-align: top;
    text-align: left;
}

#map {
    width: 100%;
    height: 100%;
}

.ol-popup {
    display: none;
    position: absolute;
    background-color: white;
    -moz-box-shadow: 0 1px 4px rgba(0,0,0,0.2);
    -webkit-filter: drop-shadow(0 1px 4px rgba(0,0,0,0.2));
    filter: drop-shadow(0 1px 4px rgba(0,0,0,0.2));
    padding: 15px;
    border-radius: 10px;
    border: 1px solid #cccccc;
    bottom: 12px;
    left: -50px;
    height: 100;
    width: 200;
    min-width: 100px;
    max-height:400px;
    overflow-y: scroll;
}

.ol-popup-closer {
    text-decoration: none;
    position: absolute;
    top: 2px;
    right: 8px;
}

.ol-popup-closer:after {
    content: "X";
}

.ol-attribution a {
    text-decoration: none;
    color: #666;
    font-family: calibri;
    font-size: 50%;
}

#popup-content>ul>li:nth-child(even) {
    background-color: #eee;
}

#popup-content ul {
    list-style-type: none;
    padding-left: 0;
}

#popup-content li {
    margin-bottom:0.25em;
}


texts: QGIS pgRoutingLayer database connection missing on localhost

texts: Copied the network tables to localhost. Now pgRoutingLayer will not show the connection under "Database" anymore. What am I missing? The username and password is saved in QGIS under the localhost connection, so it's not that.


texts:
I use QGIS 3.4. I would like to merge 2 rasters (.TIF) (using Raster-Miscellaneous-Merge) that contain elevation data. Unfortunately, I always get the error message:

  The following layers were not correctly generated.

I attach a screenshot of the log. I already checked that the files are in their correct place and played around with some of the options.

Unfortunately nothing worked. The files are enormous in size

Would this be the problem?


texts: Creating own ribbon with MapCad tools in Mapinfo Pro

texts:
I just updated Mapinfo Pro from 15 (32bit) to 16 (64bit). At work, I use the Mapcad tools. I tried to create my own ribbon with only Mapcad tools, but I don't have them in the Customize Ribbon dialog box. Or I cann't find them.

Does anyone know where to find them?

For example, I cann't find tool "update column with text object".


texts:
I've used jquery-ui draggable option as following but the moving popup can not be ended once dragging has started. Is there any another method to make OpenLayers.Popup draggable?

popup = new OpenLayers.Popup.AnchoredBubble(
    "JInfo", // id
    jFeature.geometry.getBounds().getCenterLonLat(), // lonlat
    new OpenLayers.Size(ISA.Size.POPUP_WIDTH,ISA.Size.POPUP_HEIGHT),
    html, // contentHTML
    null, // anchor
    true, // closeBox
    function(){
        // closeBoxCallback
        jSelectControl.unselect(jFeature);
    }
);
jFeature.popup = popup;
map.addPopup(popup);
$("#JInfo").draggable();


texts:
I would like to run a QgsTask when I push a button in a Qgis plugin (3.10).
My task is a subclass similar to this example
In the main class of my plugin, I have added this (connected to the button that works) :
def myBtnAction(self) :
    from .my_task import MyTask
    task = MyTask('run my task', 20)
    QgsApplication.taskManager().addTask(task)

The task is never executed.


texts: Merging adjacent and overlapping polygon from different sources with Postgis

texts:
I have three different sources in my PostGIS database that I would like to combine in a materialized view to make a single one. I want to merge all the adjacent and/or overlapping polygons to make a single geometry for each "cluster". I could compare the expected results as what the ArcMap 'Dissolve' tool would generate.

My problem is that although the result is a valid layer for QGIS, ST_Union() outputs nothing different from the input.

The script is supposed to:

Fetch the sources in a WITH statement
UNION the three sources
to make a single table
ST_Union() the geometries to merge them
and use string_agg() to keep records from merged geometries

The SQL script looks like this:

drop materialized view mh.mh_multi_source;
CREATE MATERIALIZED VIEW mh.mh_multi_source AS (
WITH mh_cic_es AS (
SELECT id, geom, CAST(type_mh AS VARCHAR(50)), source, CAST('Canards illimités Canada, 2003' AS VARCHAR(75)) AS "source_layer"
FROM mh.mh_cic_estrie
),
--SELECT * FROM mh_cic_es
mh_geomont AS (
SELECT id, geom, CAST(classe AS VARCHAR(50)) AS "type_mh", source_1 || '/' || source_2 as "source", CAST('Géomont, 2008' AS VARCHAR(75)) AS "source_layer"
FROM mh.mh_geomont
),
--SELECT * FROM mh_geomont
mh_sief AS (
SELECT id, geom, CAST('Inconnu' AS VARCHAR(50)) AS "type_mh", met_at_str AS "source", CAST('SIEF, 4eme inventaire' AS VARCHAR(75)) AS "source_layer"
FROM mh.mh_sief_view
),
mh_union AS (
SELECT * FROM mh_sief
UNION
SELECT * FROM mh_cic_es
UNION
SELECT * FROM mh_geomont
)

SELECT CAST(ROW_NUMBER() OVER (ORDER BY id ASC) AS int) AS id, ST_Union(geom) as geom, string_agg(source_layer, '/') AS "source_layer"
FROM mh_union
GROUP BY id
)


texts: Overpass (command line): how to view the results on a map?

texts:
I've tried to play with Overpass from command line because I've to do a query quite "expensive" that Overpass Turbo can't to do

My issue is about to download all the addresses in Italy

Here you are my query (it comes from Overpass Turbo ...)

<osm-script output="json">
 <union into="_">
  <query into="_" type="node">
   <has-kv k="addr:housenumber" modv="" v=""/>
   <bbox-query e="20.01708984375" into="_" n="47.15984001304432" s="36.27970720524017" w="5.9765625"/>
  </query>
 </union>
 <print e="" from="_" geometry="skeleton" limit="" mode="body" n="" order="id" s="" w=""/>
 <recurse from="_" into="_" type="down"/>
 <print e="" from="_" geometry="skeleton" limit="" mode="skeleton" n="" order="quadtile" s="" w=""/>
</osm-script>

I've saved it in a file "CiviciItalia.txt", and then I've used the following command

wget --post-file=CiviciItalia.txt http://overpass-api.de/api/interpreter --output-document=CiviciItalia.json

All works fine and I've obtained a .json file but ...... how may I visualize it on some GIS desktop tools (es. QGIS)?


texts:
I am not so seasoned at GeoTools libs, I am trying to transform a large (near 60K features) GML file. I want to change its CRS and features coordinates accordingly. The code I have is basically this:

public void parse(String crs, InputStream collection) throws Exception {
  Configuration config = new GMLConfiguration();
  StreamingParser parser = new StreamingParser(config, collection, SimpleFeature.class);
  SimpleFeature feature = (SimpleFeature) parser.parse();
  SimpleFeatureType featureType = feature.getFeatureType();
  SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
  builder.setName(featureType.getName());
  builder.setNamespaceURI(featureType.getName().getNamespaceURI());
  builder.setCRS(CRS.getAuthorityFactory(true).
    createCoordinateReferenceSystem("EPSG:4326"));
  builder.setDefaultGeometry(feature.getDefaultGeometryProperty().
    getDescriptor().getName().getLocalPart());
  Collection<Property> props = feature.getProperties();
  for (Property prop : props) {
    builder.add(prop.getName().getLocalPart(), prop.getType().getBinding());
  }
  SimpleFeatureType buildedFeatureType = builder.buildFeatureType();
  FeatureCollection<SimpleFeatureType, SimpleFeature> outCollection = FeatureCollections.newCollection();
  do {
    SimpleFeatureBuilder build = new SimpleFeatureBuilder(buildedFeatureType);
    for (AttributeType attributeType : buildedFeatureType.getTypes()) {
      Name attrName = attributeType.getName();
      Object value = feature.getProperty(attrName).getValue();
      build.add(value);
    }
    String id = feature.getIdentifier().getID();
    SimpleFeature sf = build.buildFeature(id);
    outCollection.add(sf);
  } while ((feature = (SimpleFeature) parser.parse()) != null);
}

it takes a lot of time to complete, about a minute. Is there a better solution to deal with these big feature collections? Oh, I am bound to Java 8 SE and GeoTools version 2.7.0, migrate to newer GeoTools is not an option.

UPDATE
I am dealing with a GML not shapefile. Shapefile was transformed to GML in a previous step.


texts: Defining ITRF 92 and ITRF 2008 using Proj4 parameters in QGIS?

texts:
I know how to define a CRS, but this particular one is giving me a headache. I was unable to find the EPSG code at spatial reference.

The main reason for asking this is that raster and vector data available from INEGI (Mexico) comes in ITRF 92. Also new cartography from the same source comes in ITRF 2008.

When I load a raster or a vector layer from INEGI their CRS properties only said "unknown datum based on GRS80 ellipsoid".

INEGI says that ITRF and WGS84 are "equivalent" to make maps in scales smaller than 1:50 000. Nevertheless I usually make maps in scales greater than 1:50 000.


texts: How to download GeoTIFFs using the REST API in GeoServer 2.2?

texts:
I am using GeoServer (Version 2.2 Git Revision f5b5c35076b52d02eb9cca3fa3232bc17b5f6d80 Build Date 19-Sep-2012 18:33 GeoTools Version 8.2 (rev 704570474295e339c08d1ca140d884f23a8a03a3)) and I have a workspace named myTestWorkspace and stored with name test and type GeoTIFF.

This store has URL connection parameters as follows: file:coverages/test.tif. I would like to download *.tif file by GeoServer REST API.

I found in documentation that it is possible using URL

/workspaces/<ws>/coveragestores/<cs>/file[.<extension>]

So in my case it will be /workspaces/myTestWorkspace/coveragestores/test/file.geotiff (am I right?). But when I paste this URL I've got HTTP status 404. In documentation above is information that in this case (getting 404) are two options:

  Exceptions:

  GET for a data store that does not exist -> 404 GET for a data store
  that is not file based -> 404

My data store exists (if for example I paste /workspaces/myTestWorkspace/coveragestores/test/file.tif then I've got Unsupported format: tif in response - this is ok, because tif is not valid format, so it means that GeoServer can see my store).

But is my data store file based? I thought that if data store has type GeoTIFF, then he is file based - am I right?
Any ideas why I cannot download TIFF file?

Thanks


texts: Web Services that Converts Lat Lon to Zip Code

texts: Does anyone know of a Web Services  that can convert a Lat Lon to a Zip Code?


texts: 2D polygons with height attribute to 3D features for ArcGIS Online

texts:
I have a 2D shapefile made of many polygons representing buildings. Each polygon has a height attribute. The goal is to get these into ArcGIS Online's "Scene" feature, showing extruded buildings, sitting on one of the default 3D base map terrains.

Things I have tried that were very close to working properly:

Method 1:

Inserting the 2D shapefile into ArcScene, extruding by height in the Properties, using the 3D Analyst "Layer 3D to Feature Class" to save as a shapefile, zipping up the file, importing into ArcGIS Online, adding to my ArcGIS Online scene, and configuring the layer to sit on the terrain.

Result:

Buildings appear to be triangles with a height value. I tried the same method again to see if it was just a fluke, and the triangles were in a different configuration the second time around.

Method 2:

Inserting the 2D shapefile into ArcMap, using the 3D Analyst's "Feature to 3D by Attribute" tool, selecting the height attribute as the height field, it saves as a shapefile, zipping up the file, adding to my ArcGIS Online content, adding the content to my scene, and configuring the layer to sit on the terrain.

Result:

Much better, the buildings look right - they appear to have accurate roofs and footprints - but there are no walls. I tried again with the same method in ArcScene and it did the same thing - roofs and footprints, but no walls.

SO...

How do I turn my 2D polygons with a height attribute into 3D features (you know, boxes) that I can put into an ArcGIS Online scene?

I can't use the simple extrusion feature in ArcGIS Online's extrusion property in the "Configure Layer" section of the scene, because all of the polygons are different heights.

I have access to ArcMap, ArcScene, ArcGlobe, City Engine, ArcGIS Online, and probably lots of other things. We have 3D Analyst and most other extensions.

I don't understand coding much at all so can you give me steps to accomplish this, or things to try?


texts: Unable to execute a python script that was exported from ModelBuilder

texts:
I have a model that copies a single F.C form SDE geodatabase to a scratch geodatabase then copies it back with a new name to the same SDE geodatabase again. And this works perfectly form the model.
Now I have exported the model to python script, but unfortunately it doesn’t work and I didn’t figure out what is the problem.

Here is the Result python code

    # Import arcpy module
import arcpy

# Local variables:
Points_DBO_Communities = "Database Connections\\1.sde\\Points.DBO.Communities"
Buildings_POINTS_Select = "C:\\Users\\ahmads\\Documents\\ArcGIS\\Default.gdb\\Buildings_POINTS_Select"
Points_DBO_Buildings_POINTS1 = "Database Connections\\1.sde\\Points.DBO.Buildings_POINTS1"

# Process: Select
arcpy.Select_analysis(Points_DBO_Communities, Buildings_POINTS_Select, "")

# Process: Copy Features
arcpy.CopyFeatures_management(Buildings_POINTS_Select, Points_DBO_Buildings_POINTS1, "", "0", "0", "0")

Now Trying to run the code and I got this massage :

"arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Table: Dataset C:\Users\jamaln\Documents\ArcGIS\Default.gdb\
UrbanMasterPlans_Objection_S2 does not exist or is not supported
Failed to execute (AddField)"

What might be the issue here ?


texts: Mapillary JS: mapillary-js storytelling [RouteComponent] hangs

texts:
I'm trying to use Mapillary JS (rif. http://blog.mapillary.com/update/2016/02/10/mapillary-js.html), to build a storytelling following this sample http://bl.ocks.org/knikel/630c2d6fa37a8a0e082a.

This is my code ...

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width">
  <title>JS Bin</title>

  <link rel="stylesheet" href="https://npmcdn.com/mapillary-js@1.0.1/dist/mapillary-js.min.css"/>

<button onclick="play()">Play</button><button onclick="stop()">Stop</button>

When I try to execute it nothing is happening and no errors are in the web browser console.

Any suggestions?


texts: Why would spatial join return correct distances but wrong attributes

texts:
I am doing spatial joins on points and street (line) features to identify the closest street to each point. However, despite getting appropriate distances returned for each points, closest road, the attributes of the identified road, in some cases an incorrect road and it seems to be one road specifically that is favoured. Projections of both layers and the project are consistent and nothing is selected in either layer. Am I missing anything?

Attached image shows features joined with same identified attribute for street, but distances, which are more appropriate to the points highlighted.

And the parameters for join...

Just to clear this up - the dataset I was using had been modified by others adding in extra features, which seems to have introduced my problem. I'm not sure why or how exactly, but going back to the original set with no extra features returns the results I would expect.


texts: Monitoring Geoserver for Production Environment

texts:
Is there any opensource server monitoring is available for Geoserver ? Like watchdog and nagios?

Suppose if geoserver is stopped working whether mapmeter will automatically start it..Like the following image it should happen. Is it possible ?


texts: How to create a 3D shapefile with QGIS from values stored in the attribute table

texts:
I have a shapefile containing buildings represented as polygons (basically the area of the buildings).

In the attribute table there is a field where the elevation of every building is stored (every building/polygon as a different elevation). Elevations are stored as decimal numbers.

What I am trying to do, is to convert this 2D shapefile in a 3D shapefile using the elevations stored in the "z" field of the shapefile.

With ArcGis it is possible with the tool Feature to 3D by attribute

I cannot find anything similar in QGIS.

I tried "save as" and then include "z dimension" but it does not prompt to choose a field where elevations are stored.

I also tried the v.extrude Grass command, but none of the field in the attribute table are shown to select from as height field.

Is it possible with QGIS/Grass to do this?

My final goal is to obtain buildings as 3D objects.


texts: Is it possible to publish a map service from an SD File purely over REST (no arcpy)?

texts:
I want to automate the publishing of map services, but without using the arcpy library.

Is it possible to publish a map service from an .SD File purely over REST (no arcpy)?


texts: Forcing driving direction using edited OSM shapefile and pgRouting

texts:
I need to edit a shapefile from OSM in QGIS. I want to use pgRouting to find a route, where it uses a forced driving direction. The algorithm may use the same edge however to calculate the route further. So, it has to use the forced driving direction, but it may use the same edge in the opposite direction, (unless an edge is really one-way).

Do I just add another column, and if so, how do I go about that?

More details about the use case:

Working with waste collection. Meaning the bins are on one side of the road. The truck has a loading arm on the right side. Meaning we have to drive through the street in the right way to load. After that we may use the road to transit to another part... So, in order to load the bins the road has to temporary one way but after that may be used in the opposite direction. From historical data we know where the bins are.


texts: GeoDjango - How to secure tile map services?

texts:
I use geodjango to create and serve tiles that I usually display into OpenLayers as openLayers.Layer.TMS

I am worried that anybody could grab the tms URL and plug it into their own map without asking permission, and then consume a lot of the server's CPU and violate private data ownership. On the other hand, I want the tile service to be publicly available without login, but from my website only.

Am I right to think that such violation is possible? If yes, what would be the way to be protected from it? Is it possible to hide the url?


texts: What is the difference between "Random Forest" and "Forest-based Classification and Regression" in ArcGIS Pro?

texts: ArcGIS Pro (2.6.0) has tools to train Random Forest (named Random Trees in ArcGIS) and Support Vector Machine.  Afterwards, the tool named "Classify Raster" contains the algorithms to apply your trained algorithm to imagery.  You can choose between Random Trees and SVM.  There is a unique geoprocessing tool named "Forest-based Classification and Regression" that is supposedly an adaptation of Random Trees.  What is the difference?


texts: Preserve input attribute schema option NOT working when using Feature to Polygon in ArcGIS

texts:
I have successfully converted a polyline fc with 1,780 records into a polygon fc in ArcGIS 10.1, by using the tool: Feature To Polygon, but the output does not keep the attributes of the input fc (I have checked the optional button to do so)

Then, I found Transfering attributes from polyline to polygon features using ArcGIS Desktop?, but is this the only option? because I don't want to run two tools.

What am I missing?


texts:
I am having trouble with a field calculator statement! I am trying to reclassify an attribute column based on their values. So if value =0 then make it 14, if value < 3 make it 11, if value >3 and < 9 then make it 12 and if none of the above make it 13.
I wrote this but 12 doesnt ever show up when i know it must!
My Code:
CODE BLOCK
def reclass(x):
     if x == 0:
          x = 14
     elif x <= 3:
          x = 11
     elif (x < 4 and x > 9):
          x = 12
     else:
          x = 13
     return x
field= reclass(attribute)


texts: Connecting to Cloud SQL database in Colab runtime

texts:
I am very new to GCS and Cloud SQL, so forgive me if I am not asking this question well.
A colleague and I have been working on a machine learning workflow that generates polygons from some predictions made by our model. Now, we need to store information about these polygons in a Postgres database (hosted in Cloud SQL) to be able to query information over several years of data. The issue I'm currently faced with is actually connecting to the database. Preferably, I'd like to connect to it in the Colab environment (all of our code up to this point in the workflow has been written/run in Colab) so that I can seamlessly upload the polygon information to the db using Python. Up to this point, however, I have hit a wall with every solution I've attempted.
Is there a straightforward way that I can connect to my Cloud SQL db in the Colab environment? If not, what are some alternative approaches that I might not be thinking of at the moment (preferably using Python)?


texts: How to access/edit QGIS3 Plugin folder in Flatpak version from external software in Linux Mint 19.1?

texts:
I am glad to share my first question in this forum with you. I hope I can formulate the problem correctly.

I have been working in a QGIS3 Python plugin in a Windows system with QGIS 3.8. From today on I want to develop the plugin in a computer with QGIS 3.6.2 and Linux Mint 19.1 as a OS.

I copy-pasted the plugin folder into the corresponding user profile directory of the Linux computer ("/home/USER/.var/app/org.qgis.qgis/data/QGIS/QGIS3/profiles/default/python/plugins/PLUGIN"). The plugin is working perfectly.

I could copy-paste the plugin folder and I can edit its files from the file explorer or text editors. However, I cannot access or edit the plugin folder and files from other software, specifically from PyCharm and GitKracken.
I realized that there are three subfolders in "/home/USER/.var/app/". These are: "com.axosoft.GitKraken", "com.jetbrains.PyCharm-Community" and "org.qgis.qgis". The last one is where the plugin is located.
When browsing the system from these two programs, I can only see their respective folders in "/home/USER/.var/app/", but no sign of the existing "org.qgis.qgis", as if read permissions were missing. I am working as admin user.

I would like to continue developing the plugin code with both software while it is still in the QGIS user profile folder, in order to see and reload the changes dynamically.

Any idea how to solve this?


texts: Export Polygons as Polygons and not MultiPolygons in shapefile

texts:
I have a problem when trying to export my Polygons with geopandas in shapefiles as Polygons. When I import the file later on in QGIS, they are analised as MultiPolygons.

I believe it is a tendency with the ESRI shapefile format, but is their a way around that permits to keep on with the shapefile format ? Or do I need to use another exporting format?

Here are some screenshots :

Exporting with geopandasin Python. There appears to be no problem.

Importing in QGIS. The features are imported as MultiPolygons instead of Polygons.

EDIT

The objective is to import this data to a PostGIS database using shp2pgsql. QGIS is a way to ensure that the shapefile indeed contains Polygons and not MultiPolygons. See screenshot below :

Testing if there are Polygons in the shapefile :


texts: Raster Calculator Iteration in ModelBuilder?

texts:
I have a large number of rasters and I want to add 1 to each of them (e.g. each raster's values currently range from 0 to 1, and I want to make them range from 1 to 2). I believe ModelBuilder, using the "Iterate Rasters" tool and the "Raster Calculator" tool, should do this.

In "Raster Calculator" outside of modelbuilder, selecting the raster name and typing + 1 works to do this on a single raster. However, in modelbuilder, using ("%Raster%") + 1 does not work (see attached photo).

I'm new to ModelBuilder and Python.


texts: PyQt QFileDialog custom proxy filter not working

texts:
This working code brings up a QFileDialog prompting the user to select a .csv file:

def load(self,fileName=None):
        if not fileName:
            fileName=fileDialog.getOpenFileName(caption="Load Existing Radio Log",filter="csv (*.csv)")[0]
  ...
  ...

Now, I'd like to change that filter to be more selective.  The program saves each project as a set of three .csv files (project.csv, project_fleetsync.csv, project_clueLog.csv) but I only want the file dialog to display the first one (project.csv) in order to avoid presenting the user with too many choices when only a third of them can be handled by the rest of the load() function.

According to this post, it looks like the solution is to use a proxy model.  So, I changed the code to the following (all of the commented lines in load() are things I've tried in various combinations):

    def load(self,fileName=None):
        if not fileName:
            fileDialog=QFileDialog()
            fileDialog.setProxyModel(CSVFileSortFilterProxyModel(self))
#           fileDialog.setNameFilter("CSV (*.csv)")
#           fileDialog.setOption(QFileDialog.DontUseNativeDialog)
#           fileName=fileDialog.getOpenFileName(caption="Load Existing Radio Log",filter="csv (*.csv)")[0]
#           fileName=fileDialog.getOpenFileName(caption="Load Existing Radio Log")[0]
#           fileDialog.exec_()

...
...

# code for CSVFileSortFilterProxyModel partially taken from
#  https://github.com/ZhuangLab/storm-control/blob/master/steve/qtRegexFileDialog.py
class CSVFileSortFilterProxyModel(QSortFilterProxyModel):
    def __init__(self,parent=None):
        print("initializing CSVFileSortFilterProxyModel")
        super(CSVFileSortFilterProxyModel,self).__init__(parent)

    # filterAcceptsRow - return True if row should be included in the model, False otherwise
    #
    # do not list files named *_fleetsync.csv or *_clueLog.csv
    #  do a case-insensitive comparison just in case
    def filterAcceptsRow(self,source_row,source_parent):
        print("CSV filterAcceptsRow called")
        source_model=self.sourceModel()
        index0=source_model.index(source_row,0,source_parent)
        # Always show directories
        if source_model.isDir(index0):
            return True
        # filter files
        filename=source_model.fileName(index0)
#       filename=self.sourceModel().index(row,0,parent).data().lower()
        print("testing lowercased filename:"+filename)
        if filename.count("_fleetsync.csv")+filename.count("_clueLog.csv")==0:
            return True
        else:
            return False

When I call the load() function, I do get the "initializing CSVFileSortFilterProxyModel" output, but apparently filterAcceptsRow is not getting called: there is no "CSV filterAcceptsRow called" output, and, the _fleetsync.csv and _clueLog.csv files are still listed in the dialog.  Clearly I'm doing something wrong...?


texts: I want to add a basemap under my 3D bar chart of density. I tried a lot by couldn't add the basemap in ArcScene. I want a map like the one below.


texts: Convert DateTime field to Just Date in Label Expression

texts: I have a DateTime field (ex: CreatedDate) and would like to only show the date with the label.  How would I setup the label expression to do this using python?


texts: Does openJump or Jump support the embedding of videos or images?

texts:
I've been tasking with writing a plugin for openJump that does precisely that. The idea is that you have some points/markers and when you click on them a video/image is presented (whether in a dialog, web browser, etc.)

I've been looking for other plugins, or even info on the subject for nearly two weeks and I'm beginning to think that is not possible.

Can anyone kindly provide further information?


texts: Filling table matches multilayer from buffer in PostGIS?

texts:
I have some tables in Postgis of differents geometries (Points, Lines, Polygons).

I create a buffer of 1 kilometer of one point with this:

ST_Buffer(ST_GeomFromText('POINT(3056563.45195278 804692.195990143)',48402), 1000) AS buffer

What i want to do is populate a new table with the results of my layers(tables in postgis) that are inside of the buffer.

I test with this query and it works, but only if found matches in all the tables.

    SELECT * FROM
         table1,
         table2,
  ST_Buffer(ST_GeomFromText('POINT(3056563.45195278804692.195990143)',48402),1000) AS buffer
    WHERE
      ST_Within (buffer, table1.geom) AND
      ST_Within (buffer, table2.geom)

The idea is give some point or polygon with a buffer in the map. And return a table with all the matches from table1, table2, table3, table4 that are inside that buffer.

How can i solve that problem?


texts: Plotting realtime .nmea files with QGIS 2.18.10

texts:
I am working on a project that requires us to track 2 GPS targets, one car and one weather balloon. The realtime GPS tracking of the car via the COM port works great. The decoded data from our weather balloon however is poured into a growing .nmea file, is it possible to make QGIS (2.18.10) continuously update this imported vector layer in order to display this position realtime?

So far when we import the .nmea file, only the positions up to that point are shown and the vector layer isnt updated despite the sourcefile growing in size.


texts: How to share *.sql files between ArcGIS and MapInfo?

texts:
In Denmark many municipalities have the same provider of public data. These data are pushed to a local database – often a SQL Server.
The DBMS layout is the same for all municipalities.

When I create a linked table from the SQL Server, I can save (or load) a .sql file in the MapInfo dialog.

This .sql-file is portable and I can share it with municipalities with the same setup – but only with MapInfo-municipalities.

Recently I had a post at a conference on this item and I urged my colleagues to use my .sql-files and/or my source code. I was however unable to guide my colleagues from the ArcGis-municipalities.

Can anyone give me some heads up on this matter? Is it possible to do In ArcGis?
The .sql file from MapInfo contains the plain SQL query which is sent to the SQL Server:

I.e. :

Select
"JY72000V"."VEJ_NAVN"+ ' '+ "JY72000V"."HUS_NR" As Adresse,
"JY70100V"."KOMMUNENUMMER" AS Kommunenummer,
"JY70100V"."EJENDOMSNR" AS Ejendomsnummer,
Cast("JY72000V"."KOMMUNENUMMER" As varchar(3)) + Cast("JY72000V"."VEJKODE" as varchar(4)) + "JY72000V"."HUS_NR" As ADRESSE_ID
From "ql5020pd"."dbo"."JY70100V", "ql5020pd"."dbo"."JY72000V"
Where "JY70100V"."ADG_ID"="JY72000V"."ADG_ID"

Thanks for any help

Søren Breddam,

Stevns Kommune


texts: Changing to alternate label if first label does not fit in QGIS?

texts:
I'm trying to create a map to print out and the size of my polygons varies greatly. As this will be a reference map I need to label every polygon, but some of the names will not fit in the polygon. In those cases I'd like to use the feature id instead and then have an attribute table link the id to the name.

What I would like to do is have some sort of expression that I can use to select which label to use on the map and which rows to include in the attribute table. Something like isLabelled that would be able to take the map scale, label length, font size etc into account.


texts: OpenLayers 3: Sync/Unsync Different Maps Events

texts:
Let's say I have N different maps on the screen, I want to be able based on a condition, let s say a value coming from a checkbox or a different control to be able to sync all the map events (zoom/pan/rotation) and unsync them accordingly.
 Can this be achieved in OpenLayers 3 using native API?


texts:
Most of my layers in Qgis are shapefiles but I have some raster data which is terrain relief shadings that are GeoTIFFs.

In order to save a lot of time with the default colour and density that Qgis sets, I am trying to find out if there is some way I can set the colour and transparency all at once.

Simply put I want a colour (or rendertype) of singleband grey and a global transparency of 80% to be able to be applied all at once. I tried to set up a style for this but some of this information is not saved into the style.


texts:
Working on an .net arcobjects app that's connecting to a Geodatabase (Oracle) through ArcSDE. I make a bunch of edits inside a EditSession but the data I am touching have a bunch of foreign Key constraints that I have to enforce by inserting/updating a bunch of records (non spatial) before I commit/flush/close my edit session.

Is there a way to execute those queries through the same session being used by ArcSDE without having to deal with Table feature classes? Something that allow me to run pure SQL through an edit session?


texts: Point in polygon analysis fails with "TypeError: arguments did not match any overloaded call"

texts:
I am using QGIS 1.8 on Windows 7

I want to perform a point in polygon analysis using the Vector > Analysis Tools menu.

My point layer is the result of a SpatiaLite query which returns 10'000 points.

My polygon layer is a simple shapefile with about 150 different polygons. Both layers are in WGS84 (EPSG:4326).

When I click the "OK" button in the Point in Polygon tool, just nothing happens. The progress bar stays at 0%, but I cannot click "Cancel" (i.e. nothing happens). I can click on the Close-Button in the upper right though, the window then disappears.

EDIT: This is the output from the message log:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:\OSGeo4W\apps\qgis\python\plugins\fTools\tools\doPointsInPolygon.py", line 215, in run
    tmpGeom = QgsGeometry(pntFeat.geometry())
TypeError: arguments did not match any overloaded call:
  QgsGeometry(): too many arguments
  QgsGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType'

Python version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS version:
1.8.0-Lisboa Lisboa, f440b60

Python path: [...]

I have also tried saving the point layer coming from a SpatiaLite DB as a shapefile and reimporting this shapefile. I have also tried to do it with a much smaller subset of points. No success.


texts:
After hours of tries and errors I can't get this right in QGIS 2.99.

I want to add a layer with the delimitedtext function in QgsVectorLayer but if I
add this code to my plugin, or if I'm running it in the python-console:

import os
from qgis.core import QgsProject, QgsVectorLayer

fullname = "C:/åäö/myfile.txt"
filename = os.path.splitext(os.path.basename(fullname))[0]
uri = 'file:///%s?crs=%s&delimiter=%s&xField=%s&yField=%s&decimal=%s' % (fullname, 'EPSG:3006', ',', 'East', 'North', '.')
layer = QgsVectorLayer(uri, filename, 'delimitedtext')
if layer:
    QgsProject.instance().addMapLayer(layer)
    print("Done")
else:
    print("Error")

I'm getting an error that the file can't be opened, if I change the file path to "C:/abc/myfile.txt" everything's working fine. In QGis 2.18 i solved it with adding a .encode('utf-8') to the path, but that's not working anymore.

If I'm saving that code in a file and makeing my plugin run that file with this code, it's working.
But in this way I can't use it because the file is coming from a user input through the filedialog window.

with open(u'C:/åäö/TESTdelimitedText.py') as f:
    exec(f.read())
return

or if I run it in the editor in python-console in Qgis.

In the ogr provider its working without any problems.

How do I make QgsVectorLayer delimitedtext accept a file path with special characters?


texts: Plotting raster cropped to polygon edges in R

texts:
I am trying to produce a map (using R) that uses a color-scale to display raster data.  I want to display the raster data only underneath a particular polygon.  My first thought was to use raster::mask to generate a new raster that is cropped to the polygon.  The problem is that, because of the low resolution of the raster, I end up with jagged rather than smooth edges.  Here's an example:

library(sp)
library(raster)

# Create Raster
r <- raster(ncol=40, nrow=20)
r[] <- rnorm(n=ncell(r))

# Create Polygon
a <- seq(-150,150,1)
b <- 50*a^2/(150^2)
df <- data.frame(a=a,b=b)
srl <- Polygon(df, hole=as.logical(NA))
Srl <- Polygons(list(srl), 1)
poly <- SpatialPolygons(list(Srl), pO = 1:length(Srl), proj4string=CRS(as.character(NA)))

# Plot raster and polygon
plot(r)
plot(poly, add=T)

# Mask and plot
rr <- mask(r, poly)
plot(rr)

What I want is something that looks like this instead (this code is clumsy, but it shows what I want; basically it works by plotting a big white polygon on top of the raster with a hole over the area I want visible):

# The previous code snippet must be executed prior to this one
e1 <- c(-180,180,180,-180)
e2 <- c(90,90,-90,-90)
df2 <- data.frame(e1=e1,e2=e2)
srl2 <- Polygon(df2, hole=F)
Srl2 <- Polygons(list(srl2), 2)
poly2 <- SpatialPolygons(list(Srl2))
srl3 <- Polygon(df)
Srl3 <- Polygons(list(srl3), 1)
poly3 <- SpatialPolygons(list(Srl2,Srl3))

hole <- poly
polyX <- poly2
coordsHole <-  hole@polygons[[1]]@Polygons[[1]]@coords
newHole <- Polygon(coordsHole,hole=TRUE)
listPol <- polyX@polygons[[1]]@Polygons
listPol[[length(listPol)+1]] <- newHole
punch <- Polygons(listPol,polyX@polygons[[1]]@ID)
new <- SpatialPolygons(list(punch),proj4string=poly@proj4string)

plot(r)
plot(new, col="white", border="NA", add=T)

In my case, I can't get away with this dumb hack, because the raster needs to be plotted on top of an existing plot.  So how can I plot the raster, where pixels along the polygon edge are partially plotted, so only the parts inside the polygon show?

I don't want to just increase the raster resolution with raster::disaggregate for a number of reasons.  The raster is already big to begin with, and it could get pretty unwieldy.  Also, for visualization purposes I want to interpolate the plotted raster (i.e. using raster::plot(r, interpolate='bilinear'), and I'm pretty sure disaggregate would mess with that.

Note that if I use ggplot2 to plot the raster with geom_raster, it handles the edges exactly how I want.  Unfortunately, I need to plot the raster as part of a larger plot for which I cannot use ggplot2.


texts: Trouble compiling ArcObjects Server Object Extension project into .SOE file

texts:
I'm writing a Server Object Extension(SOE) in Visual Studio 2012 using C# in Windows 7. I'm using the REST Template that comes packaged with ArcObjects 10.2.

I'm able to build the program without any errors, but the .SOE file isn't created, only the .dll file. Does anyone have any ideas on how to create the .SOE file?