texts: ArcGIS queryTask() always returns an empty set

texts:
Using the ArcGIS JavaScript API 4.14 to query a feature service on a server that is federated to my instance of Portal/ArcGIS Enterprise.  No matter what I do the queryTask.execute() is returning an empty array.

I'm using IWA to authenticate to my Portal and OAuth2 to generate a token, which queryTask() is appending to the request which looks like this:
<URL>/query/f=json&returnGeometry=false&spatialRel=esriSpatialRelIntersects&where=1%3D1&token=[tokengibberish] (that's an example where I was desperate and using a 1=1 in the where clause, just trying to get it to return SOMETHING).

My script uses hitTest() to pull an OBJECTID from a feature on the map and stores that OBJECTID to a variable:

view.on("click", function (event) {
    view.hitTest(event).then(function (response) {
        if (response.results.length) {
            let graphic = response.results.filter(function (result) {
                return result.graphic.layer === station;
                    })[0].graphic;

            // do something with the result graphic
            console.log(graphic.attributes.OBJECTID);

            document.getElementById("displaytext").innerHTML = graphic.attributes.OBJECTID;

            feature_objectid = graphic.attributes.OBJECTID;
        }
    });
});

Then I have a function built to listen for a button click event that executes the queryTask:

view.when(function() {
    document.getElementById("querybutton").addEventListener("click",
        function(event) {

            var queryTask = new QueryTask({
                url: station.url  // Station is the FeatureLayer that is displayed on the map and the OBJECTID is pulled from using hitTest()
            });

            var query = new Query();
            query.where = "OBJECTID=" + feature_objectid;  //feature_objectid comes from the result of the hitTest()

            console.log("Query: " + query.where);

            queryTask.execute(query).then(function(results){
                //Down here I'm just trying to understand the output
                console.log('Query',query);
                console.log(results);
                console.log(results.features);
            });
        }); // Close event listener
 }); // Close button view.when

No matter what I put into query.where the results are always {"layers":[]}

The URL that is sent to the server looks like this:
https://<serviceURL>/FeatureServer/query?f=json&returnGeometry=false&spatialRel=esriSpatialRelIntersects&where=OBJECTID%3D64201&token=[tokengibberish]

I can successfully run a query from the /FeatureServer/query/ page by writing in the Layer Definition manually.  I have even had success using PostMan to generate a successful query by pointing at the REST endpoint and manually copying in the token and setting the layer definition - so I know things are working I just can't seem to get my code in order.


texts: Using max cell value as input in Raster Calculator of ArcGIS Pro

texts:
I am using  ArcGIS Pro 2.4.

In the following equation, can "X" be replaced with the maximum cell value of another raster?  Alternately, replacing "X" with a field from a table would work.

Con("udlen", 1, 0, "Value > X")


texts: Overwrite the mosaic raster generated using gdal_merge.py?

texts:
Is it possible to overwrite the raster generated using gdal_merge.py? When I ran the command:

% gdal_merge.py -init 255 -o out.tif in1.tif in2.tif

it works fine, but if I ran the same command and out.tif was already previously generated, it overwrites out.tif to a uniform raster (i.e., all values equal to 255).


texts: Converting MODIS LAI data (2001) from HDF to GeoTIFF?

texts:
I am trying to convert MODIS Leaf Area Index (LAI) data from .hdf to GeoTiff from 2001.

I can successfully convert it using R for the more recent datasets using the following code:

library(gdalUtils)
sds <- get_subdatasets("GLASS01B01.V03.A2014001.2016025.hdf")
gdal_translate(sds[1], dst_dataset = "GLASS01B01.V03.A2014001.2016025.tif")

Datasets from 2014 are available here: ftp://ftp.glcf.umd.edu/glcf/GLASS/LAI/MODIS/0.05D/2014/.
I am then displaying the output GeoTIFF files in ArcGIS.

However, the result is incorrect for the data from 2001, available here:
ftp://ftp.glcf.umd.edu/glcf/GLASS/LAI/MODIS/0.05D/2001/

The spatial extent information seems to be incorrect. The spatial extent is specified within the files as Top: 0.025, Left: -0.05, Right: 0.05, Bottom: -0.025

For datasets from 2014, projected correctly the spatial extent is specified as Top: 90, Left: -180, Right: 180, Bottom: -90.

I have also tried the HEGTool.

Can I correct for this, either in R or in ArcGIS? Maybe some addition to the code will allow me to specify the extent?


texts: Calculate text field using python script and keeping preceding zeroes

texts:
I'm trying to calculate records in a text field to "01", keeping the preceding zero in tact using python script in ArcGis.

I can do this easily using the field calculator in ArcMap or the calculate field gp tool. But the script results in "1" instead of "01" and I can't figure out how to fix it.

arcpy.CalculateField_management("Table1","Field1","01","PYTHON","#")


texts:
I would like to load new composers from template files (*.pqt) when I open my plugin.

I managed to do it using this helpful post
Programmatically import new composer from template to project using pyQgis

But as explain in QGIS API doc, when using createNewComposer, the composer window will be shown and activated.

Is there a way to automatically close the composer window, then loading template would be 'transparent' for plugin user ?


texts: How to connect lines with gaps to continuous line?

texts: and I have a lot of this problem like in image. The two line are separated but i want a continuous line. What rule should I choose in Topology Checker? I want to know all the separated line, so i can fix them.


texts: Why don't I see my added layers when I choose another base map?

texts:
I'm working on a small application with OpenLayers 3. In my file.js, I have added these four base maps as shown below:

var mapquest = new ol.layer.Tile({
        source: new ol.source.MapQuest({layer: 'osm'}),
        visible: false,
        name: 'mapquest'
});
var bingmaps = new ol.layer.Tile({
    preload: Infinity,
    source: new ol.source.BingMaps({
      key: 'dY07ihXoIp15X3uz2p7u~XhaOOszGD1yRh0iA2siECg~ApUkGF8SU2HBMsmFOMX6wuknO-ehbjd919j8pD6Jg-fWqeJL97x3K3KtdApb7xTH',
      imagerySet: 'AerialWithLabels',
    }),
    visible: true,
    name: 'bingmaps'
});
var esri = new ol.layer.Tile({
    source: new ol.source.XYZ({
        attributions: [
            new ol.Attribution({
                html: 'Tiles &copy; ArcGIS'
            })
        ],
        url: 'http://server.arcgisonline.com/ArcGIS/rest/services/' +
        'World_Topo_Map/MapServer/tile/{z}/{y}/{x}'
    }),
    visible: false,
    name: 'esri'
});
var stamen = new ol.layer.Group({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.Stamen({
                layer: 'watercolor'
            })
        }),
        new ol.layer.Tile({
            source: new ol.source.Stamen({
                layer: 'terrain-labels'
            })
        })
    ],
    visible: false,
    name: 'stamen'
});

I've made visible: true for bingmaps as you can see. Then, in my file.html, I have made a list of these base maps so that the user can move from one to any other base map like he wants, as shown below:

<input type="radio" name="layer" value="bingmaps" checked>Bing Maps

<input type="radio" name="layer" value="mapquest">Mapquest OSM

<input type="radio" name="layer" value="esri">ESRI Maps

<input type="radio" name="layer" value="stamen">Stamen

The Base map bingmaps is checked.
After that, I converted four shapefiles to GeoJSON files, as shown:

var vectorSourceCentre = new ol.source.Vector({
    url: 'mysql_geojson/CENTRES.geojson',
    format: new ol.format.GeoJSON(),
    projection: 'EPSG:4261'
});

var vectorLayerCentre = new ol.layer.Vector({
    source: vectorSourceCentre,
    style: centreStyle
});

It's the same thing for the other shapefiles. I added them to the map:

var map = new ol.Map({
    target: 'map',
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine(),
        new ol.control.ZoomSlider()
    ]),
    renderer: 'canvas',
    layers: [mapquest, bingmaps, esri, stamen, vectorLayerCentre, vectorLayerPba, vectorLayerHectometre, vectorLayerChambre],
    view: view
});

When I open the browser, I see my added layers on the base map bingmaps with no problem (bingmaps is what I've chosen as visible), but the problem is when I check any other base map such as ESRI maps or Stamen I don't see  my GeoJSON layers any more.

How could I fix this?


texts: Retrieving info from GRIB2 file with JAVA?

texts:
I've just started parsing some info inside GRIB2 file with NOAA weather forecast, and I am realizing this task is not as easy as I could think 'a priori'. I need to process some info inside JAVA routines, and until now I am experiencing with library from Unidata.ucar.edu.

The problem is that accessing any specific data inside the file is a little pain in the ass, because of I had never seen a weather interchange file format until now. Therefore, while I read up documentation about GRIB2 file format, if any of you know some JAVA working examples that extract concrete information it would be very appreciated (e.g. retrieve temperature at some level in specific lat/lon)?

I found this library's API just here, but no examples are provided so in the meanwhile it's time to study.

I create a Grib2Input object with a RandomAccessFile as param like this:

Grib2Input grib2Input;
String pathfile = "c:/gfs.t06z.pgrb2f03";

//Create RandomAccessFile
RandomAccessFile gribfile = null;
try {
  gribfile = new ucar.unidata.io.RandomAccessFile(pathfile2, "rw");
} catch (IOException e) {
  e.printStackTrace();
}

//Create grib2Input file
grib2Input = new Grib2Input(gribfile);
try {
  grib2Input.scan(false,  false);
} catch (IOException e) {
  e.printStackTrace();
}

And here I can get info about GRIB2 file as Product Status, Grid Name, Grid Units, and so foth by getting Grib2GridDefinitionSection of my grib2Input object and then Grib2GDSVariables with getGdsVars() getter on it .

Now I suppose (maybe i'm wrong) each record is compound ofall sections (1, 2, 3... of GRIB2 spec) First of all: Is this righg? For example, record 1 has its GridDefinitionSection, its ProductDefinitionSection, its DataRepresentationSection, and so foth.

So, then I try to get one particular record in this way:

List<Grib2Record> records = grib2Input.getRecords();
Grib2Record myRecord = records.get(15); //Get record number 15.

How can then retrieve information just about this record as Parameter Name, Parameter Units and of course, its data??

The hole file has information about Temperature, Wind, Clouds, and much more, and my problem comes in what I'm seeing that information for each record is the same (no matter how many records I get). Of course, I'm doing something wrong or I don't know how the file is format yet.

Here are my thoughts:

1) With my RandomAccessFile, I can create a Grid2Data:

Grib2Data gd = new Grib2Data(gribfile);

2) Then I can get data in this way:

float[] data = null;
for (int i = 0; i < records.size(); i++)
{
    //First get record information
    System.out.println("record " + i +  "-> Param. Categ.: " + records2.get(i).getPDS().getPdsVars().getParameterCategory());
    System.out.println("record " + i +  "-> Param. Number: " + records2.get(i).getPDS().getPdsVars().getParameterNumber());
    try
    {
        data = gd.getData(
                records.get(i).getGdsOffset(),
                records.get(i).getPdsOffset(),
                records.get(i).getId().getRefTime() //This param must be new, because it
                                                    //is not documented in API
               );
    }
    catch (Exception e)
    {
        e.printStackTrace();
    }

    //Second print some data
    for (int j = 0; j < 50000; j+=5000)
    {
        System.out.println("data[" + j + "] -> " + data[j]);
    }

    System.out.println("--");
}

Key question:

If I open GRIB2 file with some software like tkdegrib then I can see information about a number of records (pressure, temperature, etc) that not match with the number of records I can see with the routine above. For example above, with a GRIB2 from NOAA about weather forecast, tkdegrib shows information about 319 records and I can print information of 315 records ¿? So I can't match my printed data with the extracted with tkdegrib to test this information is the right one.

Any suggestion?

Am I doing something wrong?


texts: R extract function buffered points weighted mean not returning the right value

texts:
I want to do weighted area mean of raster cells included in buffered points,spatialpointsdataframe but from what I have noticed the extract function only includes raster cells having their cell centers inside the 7 meter buffer. Raster cell size is 12.5 meter. It would be great if the function could calculate a weighted area mean of all raster cells not just cells having their centers inside the buffer zone. Anyone have an idea how to solve this?
Code looking like this:

means <- extract(Raster, SpatialPointsDataframe, method='simple', buffer=7,weights=TRUE, normalizeWeights=T, small=T, fun=mean, df=T)

As an example this spatial point gets the value of one pixel with weighted mean 21


texts: Filter a QGIS report to get only selected features information

texts:
I have a layer with about 28000 points. What I'm trying to do is create a report showing just the information of the selected points, but looking at the report manager options I have not found an option to make a filter of this type.
Is it possible to make a filter like this? Maybe programatically?


texts: How to correctly plot SpatialPolygonDataFrames in ggplot and substract them in R

texts:
I have a landuse shapefile (contains data and script) that plots differently in Base and ggplot, the reason for which remains a myster after countless hours of playing with it.

I read in both shapefiles as follows, since catchment_landuse does not have an ID column, I add one:

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, raster, broom, rstudioapi)
#setwd to file location
setwd(substr(rstudioapi::getActiveDocumentContext()$path, 1, rev(gregexpr("/", rstudioapi::getActiveDocumentContext()$path)[[1]])[1]))

catchment_industry <- readOGR('.', 'catchment_industry')

catchment_landuse <- readOGR('.', 'catchment_landuse')
catchment_landuse@data$id <- 1:nrow(catchment_landuse@data)

Is there actually a difference in

catchment_landuse@data$id <- 1:nrow(catchment_landuse@data)

and

catchment_landuse$id <- 1:length(catchment_landuse$id)

?

My code for plotting the thing in base:

base_plot_spdf <- function(spdf){
  # get classes for for-loop and sort alphabetically
  classes <- sort(unique(spdf$Class))

  # rainbow color vector
  class_colors = rainbow(length(classes))

  spdf_bbox <- bbox(spdf)

  # enable legend outside of plot
  par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)

  #expanded empty plot for legend
  plot(NULL, xlim = spdf_bbox[1, ], ylim = spdf_bbox[2, ], ylab = "", xlab = "")

  # plot polygons that belong to specific classes

  for(i in 1:(length(classes))){
    plot(
      subset(spdf, Class == classes[i]),
      col = class_colors[i],
      add = TRUE
    )
  }
  # legend
  legend("topright", inset=c(-0.25,0), legend=classes, col = class_colors, lwd=c(2.5,2.5) )
}

base_plot_spdf(catchment_landuse)

gives

to make this plottable in ggplot and keep the attributes, I do

spdf_to_df <- function(spdf){
    tidy(spdf) %>%
    merge(spdf, by = 'id') %>%
    as.tibble
}

catchment_landuse_df <- spdf_to_df(catchment_landuse)

Then I use this code for actual plotting:

plot_spatial_df <- function(df){
  ggplot(df) +
    geom_polygon(aes(
      x = long,
      y = lat,
      group = group,
      col = Class,
      fill = Class)
      ) +
    scale_fill_manual(values = rainbow(length(unique(df$Class)))) +
    theme_minimal() +
    coord_equal()
  }

plot_spatial_df(catchment_landuse_df)

Which gives

as you can see, especially in the lower left and lower right classes do differ completely. What is the reason for this?
I do not need the "id" column that I create for base graphics, however for ggplot I need to do the "self-join" on the "id" column to keep the attributes.

Also, why does one id get lost somewhere?

length(unique(catchment_landuse_df$id))
[1] 575

How does the @data of the spdf know to which polygons a row in @data refer? What is the link for that?

Now I want to substract the second shapefile, catchment_industry from that, keeping all attributes in catchment_landuse:

catchment_landuse_without_industry <- erase(catchment_landuse, catchment_industry)

Base produces the output I expect:

However, ggplot now does really weird things:

catchment_landuse_without_industry_df <- spdf_to_df(catchment_landuse_without_industry)

plot_spatial_df(catchment_landuse_without_industry_df)

Is the underlying data corrupted in any way or does the mistake solely stem from wrong plotting?


texts: Propose a new location for Watch Tower using Viewshed Analysis

texts:
I want to suggest location for new watch tower in a forest area, but i am confused , on what criteria i can propose a best location. If i am taking DEM , slope and distance from existing tower , then how can i propose a best location in terms of GIS. I don't want to select points mannualy and check viewsheds , again a new point and its viewshed, until I get a maximum visible area.I want to utilize gis tools.
I tried with dem ,slope, proximity from existing tower and there reclassification into 10 classes and then weightage overlay with more percentage to dem, then distance from existing tower and least to slope.But I couldn't get desired results.
Please suggest me a solution, I will be very thankful to you friends.

Friends today I have proceeded with given problem as follows:
I want to locate a watch tower to monitor forest fire:----
1) First Criteria: Distance From Roads and Drainage ( 200 meter and 500 meter Buffer respectively), and Intersect them to find suitable areas .
2) Then I have made slope file from SRTM DEM and taken less then 6 degree as criteria (using Raster Calulator) , and I got almost 85 % study area comes under this range. DEM ranges from 400 to 763 for my study area, and using Raster Calculater I got new layer with elevation > 500 (ROUGHLY). And Converted both slope layer and DEM to polygon.
3) ROAD_DRAINAGE_BUFFER is now INTERSECTED with SLOPE_DEM Layer. And I got fewer Polygons, but Still many.
Now if I am trying to get Viewshed at highest possible point , then again viewshed from this point, is not giving me maximum visible area.
Still I am getting big area to analyse. How to get reduce my area. Still a problem.
I don't want to make you uncomfortable , asking same problem.
Please If you still give some clues, I will be thankful.
Thank You Friends.,....for giving most precious .thing......(TIME)........


texts: How to troubleshoot a hard crash on CopyFeatures_management()

texts:
I have a script that uses arcpy.CopyFeatures_management(FEATURE_CLASS_SDE, FEATURE_CLASS_LOCAL) to copy a feature class from our enterprise SDE to a local file geodatabase. For the most part, it works well. Occasionally, though, it fails in spectacular fashion.

when it fails, I get a message that states python stopped working with a minidump. That's it. I wrapped the line that performs the copy in a try...except to get more info but the exception is not caught.

Any ideas how I can troubleshoot this further?

ArcGIS 10.1 SP1 w/ QIP (Build 3300)

Update 1: When I run the script within PyCharm, I'm unable to reproduce the issue. When I call the script from a batch file, on the same workstation, the issue occurs. This workstation has 8GB RAM, plenty of HDD space available, and by watching Resource Monitor I don't see any indications we're even close to running low on resources. Python.exe is using about 300mb. There are 95k features in the this SDE FC.


texts: Converting vertical datums GeoTIFF to GeoTIFF

texts:
Anyone know where I've gone wrong here?
gdalwarp input: geotiff that is in NZVD2009 [EPSG:4440]
gdalwarp output: geotiff that is in NZGD2000 [EPSG:4167]
(gdal) C:\>gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=C:/dev/lidar-qc/las_qclegacy/vdatum_nzvd09_16/nzgeoid2009.gtx" -t_srs "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" "D:\gerald_conv\DEM_BY18_2016_1000_4646.tif" "D:\gerald_conv\DEM_BY18_2016_1000_4646_nzgd2000.tif"
ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds.
(gdal) C:\>


texts: How to plot certain values of a column in R?

texts:
Lets say we have a shapefile of cities with a column called category

cities@data

category
 A
 B
 A
 C
 A
 A
 B
 C

How can i plot the A category only etc?

Tried

  c<-cities$category[A]
  plot(c)

but wasn't correct.


texts: How to convert ASCII in x,y,z format using gdal_grid?

texts:
I have a map saved in (AFAIK) Arc/Info ASCII format and I would like to convert it with gdal grid but I am getting errors.

My data (p.csv) is in the following format:

319053.25,319301.12,544.28
319074.57,319275.26,546.44
319075.00,319250.00,546.34
319075.00,319000.00,468.74
319075.00,318950.00,466.19
319075.00,318925.00,472.06

and I have a p.vrt written by myself:

<OGRVRTDataSource>
    <OGRVRTLayer name="dem">
      <SrcDataSource>p.csv</SrcDataSource>
      <GeometryType>wkbPoint</GeometryType>
      <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

When I launch gdal_grid I get the following error :

$ gdal_grid -l dem p.vrt p.tiff
ERROR 1: Failed to find layer 'dem' on datasource 'p.csv'.
No point geometry found on layer dem, skipping.


texts: Plotting data only inside polygon in MATLAB?

texts:
We can plot data with polygon in MATLAB as a data size of 651*1200*12
S.X and S.Y shows the X and Y values of polygon.

for i = 1:12
subplot(4,3,i)
pcolor(lon1,lat1,data(:,:,i)) ;
shading interp ;
geoshow([S.Y], [S.X],'Color','black');
drawnow
axis tight
end;

This method plots the all data both in and outside polygon but I have to plot data only inside the polygon. So how we can plot the data only inside polygon and put specific colour(white/black/green etc) to area outside the polygon.

The above is the output of my way.

Yet, I want to get this kind of plot:

Only plot the data inside polygon with least gap between subplots both by row and column.


texts: QGIS Georeferencer Problem After Photoshopping Source Raster

texts:
In QGIS 3.8.2 I georeferenced a set of aerial photos from USGS Earthexplorer website and decided that I need to photoshop them to brighten them. So I photoshopped the original ungeoreferenced source images and saved them at the same resolution/size.

Then I opened one of these processed images in QGIS Georeferencer and loaded the GCP points I had previously created for the unprocessed source image. The points load above the processed image as shown in the picture below.

There may be some geodata in these images from the USGS website that gets reversed somehow when I am saving them in Photoshop?

I have check the image sizes and the processed and source images are exact same pixel dimensions.

I tried georeferencing the processed image with four points, then reloading the outputted raster back into georferencer, and loading the original set of GCP's but got same result.


texts: Creating one-to-one relationships in QGIS

texts:
Can I create one-to-one relationship using project relations in QGIS?
I am trying to create a spatial database in QGIS. I want to link my attribute data to my spatial data (shapefile). I have already done this using attribute join. Do I need to go ahead and define the relationship using project relations or does the join (which I already created) create the spatial database?


texts: Geonames Search API Result does not contain all alternate names

texts:
The entity 5124497 ("Lindenhurst") has an alternate name "Breslau":
https://www.geonames.org/5124497/lindenhurst.html
However, when using the search API the alternate name "Breslau" is not included:
http://api.geonames.org/searchJSON?username=USERNAME&style=full&name=Lindenhurst
Even though the alternate name "Breslau" is not included, the entity is included in the results for a name=Breslau search:
http://api.geonames.org/searchJSON?username=USERNAME&style=full&name=Breslau
Other examples for this API behavior are:

https://www.geonames.org/3100300/dobrowo.html, alternate name: "Post
Danzig"
https://www.geonames.org/5090215/town-of-newbury.html, alternate
name: "Danzig"
https://www.geonames.org/3099433/gulf-of-gdansk.html, alternate name: "Gulf of Danzig"

It seems that if the lang attribute of an alternate name is missing,  the alternate name is never included in the API result (although the entity can be found via an API search for the alternate name).
Is this a bug in the API or is this behavior intended? (How) is it possible to get all the alternate names?


texts:
I get a windows error when trying to run my tool again in the same arc project. It does run again once I shut down ArcMap and reopen and run the tool.  Any Ideas?

  import arcpy, os, shutil, sys

# Path to store geoprocessing files
user_profile_path = os.environ['USERPROFILE']
gdb = r"{}/AppData/Local/Temp/PolyChecker/PolyCheck.gdb".format(user_profile_path)

# GDB already exists, delete it
if arcpy.Exists(gdb):
    shutil.rmtree(r"{}/AppData/Local/Temp/PolyChecker".format(user_profile_path))

# Create fresh GDB path
os.makedirs(os.path.dirname(gdb))

# Create GDB
arcpy.CreateFileGDB_management(os.path.dirname(gdb), os.path.basename(gdb))

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

MapUnitPolys_CopyFeatures = arcpy.GetParameterAsText(1)

# Set Geoprocessing environments
MapUnitPolys = MapUnitPolys

# Validate that all Polygons have a map unit

invalid_polygon_found = False

with arcpy.da.SearchCursor(MapUnitPolys, ['SHAPE@', 'MapUnit', 'OBJECTID']) as cursor:
    for row in cursor:
        #arcpy.AddMessage(str(row[1]))
        # Does this Polygon have a map unit
        if row[1] == "" or row[1] == "<Null>" or row[1] is None or row[1] is 0:
            invalid_polygon_found = True
            arcpy.AddMessage('Polygon OBJECT ID:{} is missing map unit... exiting.'.format(row[2]))

# Invalid polygons were found, terminate
if (invalid_polygon_found):
    sys.exit(1)

Polygon_Neighbors = "{}/polytest".format(gdb)

PolygonNeighbor_TableSelect = "{}/PolygonNeighbor_TableSelect".format(gdb)

inFeatures_lyr = "Polygon_Neighbors_tv".format(gdb)

# Process: Polygon Neighbors
arcpy.PolygonNeighbors_analysis(MapUnitPolys, Polygon_Neighbors, "OBJECTID;MapUnit", "NO_AREA_OVERLAP", "BOTH_SIDES",
                                "", "METERS", "SQUARE_METERS")

Polygon_Neighbors_tv = arcpy.MakeTableView_management (Polygon_Neighbors, "Polygon_Neighbors_tv") [0]
# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(Polygon_Neighbors_tv, "NEW_SELECTION", "src_MapUnit = nbr_MapUnit")

# Process: Table Select
arcpy.TableSelect_analysis(Polygon_Neighbors, PolygonNeighbor_TableSelect, "src_MapUnit = nbr_MapUnit")

arcpy.GetCount_management(PolygonNeighbor_TableSelect)

arcpy.AddMessage(arcpy.GetMessages())

if int(arcpy.GetCount_management(PolygonNeighbor_TableSelect)[0]) > 0:
    arcpy.MakeFeatureLayer_management(MapUnitPolys, inFeatures_lyr)
else:
    print ("done")

# Process: Add Join
arcpy.AddJoin_management(inFeatures_lyr, "OBJECTID", PolygonNeighbor_TableSelect, "src_OBJECTID", "KEEP_COMMON")

# Process: Copy Features
arcpy.CopyFeatures_management(inFeatures_lyr, MapUnitPolys_CopyFeatures, "", "0", "0", "0")

#arcpy.DeleteField_management("MapUnitPolys_CopyFeatures"["IdentityConfidence", "NODE_COUNT", "OBJECTID", "Label", "Notes", "Symbol"])

# Process: Remove Join
arcpy.RemoveJoin_management(inFeatures_lyr, "")

# Execute Delete
arcpy.Delete_management(PolygonNeighbor_TableSelect)
arcpy.Delete_management(Polygon_Neighbors)

arcpy.AddMessage('All done! Check Polygons')

ERROR

  Traceback (most recent call last):
      File "C:\Users\Downloads\GeMS-ArcPro-Toolbox-master\GeMS-ArcPro-Toolbox-master\Scripts\PolyChecker2.py", line 9, in <module>
        shutil.rmtree(r"{}/AppData/Local/Temp/PolyChecker".format(user_profile_path))
      File "C:\Python27\ArcGIS10.6\Lib\shutil.py", line 261, in rmtree
        rmtree(fullname, ignore_errors, onerror)
      File "C:\Python27\ArcGIS10.6\Lib\shutil.py", line 266, in rmtree
        onerror(os.remove, fullname, sys.exc_info())
      File "C:\Python27\ArcGIS10.6\Lib\shutil.py", line 264, in rmtree
        os.remove(fullname)
    WindowsError: [Error 32] The process cannot access the file because it is being used by another process: 'C:\\Users\\/AppData/Local/Temp/PolyChecker\\PolyCheck.gdb\\a00000009.gdbtable'


texts: Computing new attribute based on changes in another attribute using ArcGIS Desktop with Python?

texts:
I am attempting to classify a set of gps time-encoded point data into behaviors based on different attributes.

I have created an attribute that is 0 for home and 1 for away based on location, and now want to number the trips away from home (a set of points 01111111111110 would be one trip because it started and ended at home). I have added the attribute field that will have the trip numbers, but don't know how to calculate the field so it is based on the home/away field.

Here is an example of the GPS data (using "*" to indicate irrelevant information and simply indexing times as 1, 2, etc.), the "Home/Away" indicator described above, and the desired trip indicator, "Trip", which I need to compute:

Time Lat Lon Home/Away Trip
   1   *   *         0    0
   2   *   *         1    1
   3   *   *         1    1
....
  12   *   *         1    1
  13   *   *         0    0
  14   *   *         0    0
  15   *   *         1    2
  16   *   *         1    2
....
  34   *   *         1    2
  35   *   *         0    0
  36   *   *         0    0
  37   *   *         1    3
....

My data set is too large to manually go through and number each trip in the attribute table, so is there any way to calculate the field based on how the home/away attribute is ordered and each "clump" of away points is designated as a trip?

These are the bare bones of what Python code might look like (I'm not experienced with code).

Expression:

trip = Reclass(!home!)

Codeblock:

def Reclass(home):
  if (home = 0):
    return 0
  elif (home = 1 and lastValue = 0):
    return _(incremental numbering?)_
  elif (home = 1 and lastValue = 1):
    return lastValue

After using matt wilkie's recommended script I've made some changes so that my first trip is number 1, my second is 2, etc etc.

Here is the code modified from matt's:

import arcpy
rows = arcpy.UpdateCursor("test2")

trip = 0
for row in rows:
    if row.home == 0:
        prev = row.home
        row.TRIP = trip
        rows.updateRow(row)

    elif row.home == 1 and prev == 0:
        trip += 1
        prev = row.home
        row.TRIP = trip
        rows.updateRow(row)
        rows.next()

    elif row.home == 1 and prev == 1:
        prev = row.home
        row.TRIP = trip
        rows.updateRow(row)
        rows.next()

    row.TRIP = trip
    rows.updateRow(row)

del row, rows

Then I just select for home = 0 and calculate my trip field back to 0. Neatly ordered trips.


texts: Remove feature class from group layer in C#

texts:
I have a large program that pulls in various group layers and feature classes. What I am trying to do is check each feature class within a specific group layer that is empty. The count part is easy and just leaves the feat class unsourced. I really want to remove it when empty. This is a brief code example:

public static void FoiUpdSym()
    {
        IMxDocument MXDoc = (IMxDocument)PublicVars.g_App.Document;
        IMap map = MXDoc.FocusMap;
        IFeatureWorkspace SDEFeatWS = null;
        IFeatureWorkspace SDEMasterWS = null;
        IFeatureWorkspace FoiUpdateWS = null;
        IGxLayer GxLayer;
        IGxFile GxFile;
        ILayer layer;
        UID uid;
        IEnumLayer enumLayer;
        IFeatureLayer featLayer = null;
        IFeatureClass featClass = null;

        // load group layer file
        GxLayer = new GxLayer();
        GxFile = (IGxFile)GxLayer;
        GxFile.Path = "\\\\gisTools\\gis\\symbology\\foiupdate\\foi-update_2016.lyr"; // strSymPath
        GxLayer = (IGxLayer)GxFile;
        layer = GxLayer.Layer;
        map.AddLayer(layer);

        uid = new UIDClass();
        uid.Value = "{40A9E885-5533-11d0-98BE-00805F7CED21}"; // IFeatureLayer
        enumLayer = map.get_Layers(uid, true);
        layer = enumLayer.Next();

        while (layer != null)
        {
            featLayer = (IFeatureLayer)layer;
            if (Regex.IsMatch(layer.Name, "cityp_adds"))
                {
                    try
                    {
                        featClass = FoiUpdateWS.OpenFeatureClass("cityp_adds");
                        IQueryFilter queryFilter = new QueryFilterClass();
                        if (featClass.FeatureCount(queryFilter) > 0)
                        {
                            featLayer.FeatureClass = featClass;
                            featLayer.Selectable = false;
                            featLayer.Visible = false;
                        }
                        else
                        {
                           map.DeleteLayer(featLayer);
                        }
                    }
                    catch { }
                }

The else statement is what isn't working.


texts: How to measure RTK accuracy of a moving vehicle?

texts:
I am currently working with tractors using a RTK enhanced automated steering system. One of the biggest problems is that even RTK can have a inaccuracy of 2cm. This can be too much for some possible field applications like sowing. Therefore I am trying to determine the inaccuracy of such a steering system in dependence to the distance of the RTK-station. My problem right now is, that I don't have a way to measure the inaccuracy of the tractor driving the track.

In my current setup I want the tractor to drive 100m straight between two given points. Thereby I want to measure the inaccuracy. This I want to repeat with several distances to the RTK-station (1km, 2km, 3km, 4km, 5km, 10km, 15km, 20km, 30km, 40km). I hope that I can get a function which tells the inaccuracy in dependence to the distance of the RTK-station.

Quenstions:

Is there a way to process the points given by the GPS device to
determine the inaccuracy?
Do you know any visual validation method? (with stereo cameras or something)
Do you have any other idea of how to measure the inaccuracy ?


texts:
I am trying to reproject an XYZ RGB point cloud file that I got from a drone mapping software. It is now in WGS84 and I want it in euref32(UTM) or just a local metric system so i can import it Into an CAD program. There is an option to remove the RGB values, but then I would lose the colors in my model.

How can i do this without losing the RGB values?

I have tried the reprojecting function in QGIS and only gets error.
This is the error I get:

ERROR 1: Ungridded dataset: At line 3, change of Y direction

Here is an sample from my XYZ RGB file:

-117.259861 32.802178 -0.511405 110 114 121
-117.259718 32.802270 4.390094 116 105 84
-117.259853 32.802182 -0.488624 114 115 118
-117.259799 32.802216 1.245612 63 73 62
-117.259711 32.802274 4.404747 115 105 83
-117.259772 32.802232 1.830184 56 81 42
-117.259749 32.802247 2.655255 77 89 55
-117.259745 32.802249 2.827860 130 124 98


texts:
Whenever I try to export a map as an PNG image I get the 'successfully export' message but there is no file generated. The error happened in Zanzibar so I changed to Noosa but I still can't export the map.

It has never happened to me before so I guess it might be a problem with my system.

Does anybody have any idea of what is happening?


texts: Where to create a web map of polygons and 50,000+ points?

texts:
I was using the qgis2web plugin to create a web map of polygons. I got this far: vehiclethefts.netne.net
Then, I had to add another layer containing 50,000 points and realized the plugin couldn't handle it. So, it crashes and you can not navigate on the map due to the large number of points.

So, I was wondering if there is a similar free service that would allow me make such a map of both polygons and points.


texts: Moving points onto lines (~neighborhood)

texts:
I do have two vector layer, of which
one is a point layer based on "events" by remote sensing
and the second one is a line layer from local research.

In my case these are earthquakes and tectonic faults, but
I guess one could simply choose "car-accidents and roads"
as a general example.

So what I'd like to do is move/copy the points onto
the closest point of the lines, as long as its within a
tolerance distance (say 1-2km or 0.0xx°), with the new
point layer (+attr moved y/n).

Any ideas ?

Linux, QGIS 1.8


texts:
One of the biggest criticisms of GIS professionals is typically that GIS is an underused technology and as such the professionals feel frustrated that they only get to use it in very limited situations such as basic map generation.

What are the top unsung business benefits that can be brought by GIS? This leads onto the next logical step of "How to sell GIS to a business".


texts:
I'm searching for a certain pathfinding algorithm which considers a certain sort of branch in the network. However, I find it difficult to define my search term as I do not know how this certain kind of edge is named. Have a look at the image, the red circled branch.

How to name a branch of which one of the two nodes is only single connected, in other words, one node is integrated in the network by only one connection?

Usually, in a network, most nodes are more or less 'intersections', meaning they connect two edges. However, in the image, you see that one edge 'sticks out', like a street ending in nowhere.

How to name this sort of edge?

Sorry, if the answer to this question may be to some of you very natural, I'm very new to this topic.
This question relates to my other question ('Yen's or Eppstein for path with intermediate destinations with 'dead ends'?').


texts:
QGIS currently does not have a way to set a default font to be used for labelling, and being lazy I'm getting tired of having to set our house font styling for every layer I create or use (see change default font of labels in QGIS for someone else equally frustrated, and the bug report for it).

I've looked at the QGIS code changes required and it's not easy, so now I'm trying for an interim fix by writing a plugin that will either apply a default qml style file or directly set it via the api.

A couple of approaches I'm contemplating are:

A legend context menu entry (easy to do, expandable to choosing from a list of house styles, but still a manual process)
Catching the layerAdded signal, checking if the new layer has a style already set, and if not setting the font (more complex but fully automatic)

Before I do this, can anyone suggest a simpler hack, or something obvious I'm missing? Something using project templates, default styles, layer definition files, or something similar?

(I'm aware this won't work when using rule-based labeling, but for most use cases it will make life easier).


texts:
Some days ago I upgraded from QGIS 2.18 to 3.6.0 - not knowing that many Python scripts won't run without modifications.

I could convert all of them but not this one:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4 import QtGui
import processing, sys, time, psycopg2
import pprint

...
creatin layers skipped (got this running)
...

################### Create Graph #####################

mapRenderer = iface.mapCanvas().mapRenderer()

layers = [missingLayer.id(),clipLayer.id()]
mapRenderer.setLayerSet(layers)

c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)

c.setPaperSize(148,90)
c.setPrintResolution(300)

# add map
x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x,y,w,h)
c.addItem(composerMap)

#output to image
dpi = c.printResolution()
dpmm = dpi / 25.4
width = int(dpmm * c.paperWidth())
height = int(dpmm * c.paperHeight())

BackgroundColor = QColor(0,0,255)

# create output image and initialize it
image = QImage(QSize(width, height), QImage.Format_ARGB32)
image.setDotsPerMeterX(dpmm * 1000)
image.setDotsPerMeterY(dpmm * 1000)
image.fill(BackgroundColor)

header = QgsComposerLabel(c)
header.setItemPosition(40,1)
header.setBackgroundEnabled(0)
newFont = QFont("Comic Sans MS", 18)
header.setFont(newFont)
header.setText("Missing Countries Map")
header.adjustSizeToText()
c.addItem(header)

header2 = QgsComposerLabel(c)
header2.setItemPosition(50,11)
header2.setBackgroundEnabled(0)
newFont2 = QFont("Comic Sans MS", 10)
header2.setFont(newFont2)
header2.setText("TS="+ts)
header2.adjustSizeToText()
c.addItem(header2)

footer = QgsComposerLabel(c)
footer.setItemPosition(124,87.5)
footer.setBackgroundEnabled(0)
newFont3 = QFont("Times Roman", 4)
footer.setFont(newFont3)
footer.setText(time.ctime())
footer.adjustSizeToText()
c.addItem(footer)

legend = QgsComposerLegend(c)
legend.model().setLayerSet(mapRenderer.layerSet())
legend.setItemPosition(5,65)
legend.setFrameEnabled(1)
newFont = QFont("Comic Sans MS", 6)
legend.setStyleFont(QgsComposerLegendStyle.Title, newFont)
legend.setStyleFont(QgsComposerLegendStyle.Subgroup, newFont)
legend.setStyleFont(QgsComposerLegendStyle.SymbolLabel, newFont)
legend.setStyleMargin(QgsComposerLegendStyle.Title,2)
legend.setTitle("Legend")
legend.setSymbolHeight(2)
legend.adjustBoxSize()
c.addItem(legend)

# render the composition
imagePainter = QPainter(image)
imagePainter.setRenderHint(QPainter.Antialiasing)
sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
targetArea = QRectF(0, 0, width, height)
c.render(imagePainter, targetArea, sourceArea)
imagePainter.end()
imageFile = "missing_boundaries"+ts+".png"
print "saving to "+imageFile
image.save(imageFile)

Yes, I do know the QqsComposer has gone but I can't find any documentation or example for QGIS3 plotting. Except the API description of course.

I don't need the whole script converted - only a hint how to open/create a png the new way.


texts:
I would like to split areas on a map of north germany between four community centers to split the responsibility for the whole map between these 4 areas.

Is there an online tool where I can draw the borders and store them so I can discuss the frontier?

I created a map with Gimp now, but that seems a bit complicated to change: https://github.com/rubo77/ffnord-community-borders


texts: Getting ArcMap layer with symbologyType = "UNIQUE_VALUES"?

texts:
I am trying to write an automated mapping script to set new categories for a layer on each iteration, however I have encountered a problem.

If I understand it correctly may layer's symbologyType should be "UNIQUE_VALUES" to be able to modify the available categories. However my dataset has the symbologyType "OTHER".

Now I understand that symbologyType is a read-only attribute and I need to use UpdateLayer to create a new layer with the desired type, but where do I get a layer with the symbologyType "UNIQUE_VALUES" to do this?

All layers I have created have type other.

I am using ArcGIS Desktop 10.6.1


texts: Change polygon color over time by number of events in polygon

texts:
In ArcGIS Desktop 10.5.1, I have a layer of XY events with a date field, and a layer of polygons. I have spatial joined the points so that each event is associated with a polygon. I want to show the polygons using graduated colors representing the number of events within it in a given time window. Ideally, this would respond dynamically to the length of the time window - e.g., I could change the aggregation level from monthly to annual with just the corresponding time slider option, and not have to mess with the data - but I'll settle for what I can get at this point.  I've approached this from a few different angles, but they all appear to be blocked:

Spatial joining the polygons onto the points doesn't work, because the geometry type is still Point and not Polygon.
Calculating Frequency on the points layer by polygon-year and joining this to the original polygon layer doesn't work, because the join just picks the first row from the frequency table for each polygon.
A relation seems to give the one-to-many polygon to polygon-time mapping I want, but it doesn't seem like you can set symbology from a relation.

I found a blog post from 2011 that does exactly what I want, but the comments are mostly people saying they're unable to reproduce it because of the enforced 1-to-1 behavior for attribute joins in ArcMap, and suggesting elaborate and unscalable workarounds. I'm using a file geodatabase, so I can't use a query table/layer.

Is it possible to achieve the behavior I want without writing a custom layer/renderer, or moving to a query layer-compatible geodatabase type?


texts: PostGIS, Generating a city "code" for every city on earth

texts:
I'd like to be able to query a PostGIS table with a lnglat and have it return a four letter.

The point [180, 90] should return AAAA, and the point [-180, -90] should return ZZZZ.

Kind of like this: https://map.what3words.com/soon.escalator.minimalists

Maybe this is just a simple fishnet grid? Is that the right way to make this?

I would like the grid to be accurate to the tenths or hundredths column of latitude. It doesn't have to be a square grid.

Cities can have multiple codes but each city should have a unique code.

edit: I realize that 26^4 might not be enough to cover the whole earth with the accuracy that I desire. So I'm thinking since most of the earth is water we could do something like this:

http://blog.cleverelephant.ca/2018/06/polygon-splitting.html

maybe this isn't possible with four letters, but six letters it should be doable if we ignore ocean (361800000 sqkm / 26 / 26 / 26 / 26 / 26 / 26) right?


texts: Euclidean distance from just one feature in a raster layer or polygon

texts: I have a raster layer ( that I converted in polygon as well) in which there are villages (black in the picture). I would like to do a euclidean distance but just from this feature "village". I was thinking to extract just this feature but I don't know how to do.


texts: Writing shapefile field values to txt file using SearchCursor?

texts:
I'm having trouble finding a way to write the output of a SearchCursor function to a txt. file. I can print it just fine in a format that is desirable. I just don't know how to get from there to a txt. file. Here's my code below:

import arcpy
from arcpy import env
env.workspace = "G:\Geocomputation_Project\Section C\Lin_Ref_test"
shape = "Points.shp"
cursor = arcpy.SearchCursor(shape)
f = open("test_field.txt", "w")
for row in cursor:
    print (row.getValue("ITEMID")),(row.getValue("New_X_LRR")),  (row.getValue("New_Y_LRR"))
print "done"

The txt file does not yet exist. When I start running this script, I hope to be able to overwrite the text file as I do more iterations.

I am using ArcMap 10.2.2


texts: Masking land with matplotlib basemap in Python?

texts:
I've been trying to mask land with .is_land from mpl_toolkits.basemap.

When running following code:

map = Basemap(llcrnrlon = 7.298914095230308, llcrnrlat = 58.98235690451632,
           urcrnrlon = 12.27072348324015, urcrnrlat =
           57.92306182768044,projection='cyl', resolution='f')
value = map.is_land(11.61168822665539, 57.86868795294363)

print(value)  returns False

but when not specifying region:

map = Basemap(projection='cyl', resolution='f')
value = map.is_land(11.61168822665539, 57.86868795294363)

print(value)  returns True

I can't figure out why.

I really do need to specify region or else the code will be 5-6 times slower to run.


texts: Why is Update Cursor adding Object IDs incorrectly?

texts:
I am using an update cursor to add values in a try:except block of code using ArcPy.  If the code is successful I log yes in the success field, if it reaches the except block I log yes in the fail field.  I also log the time of the succeeds or fail.  However, my Object ID is doing very strange things.

The first two object IDs populated this table on one execution of the script, the others came at 2 other times.  I would ideally like to see 1-6, not 1,2,401,402,801,802 as noted in the screenshot.

code snippet which reflects table

except:
rows = arcpy.InsertCursor(aTable)
for x in xrange(0, outputobjects):
    row = rows.newRow()
row.setValue("Success", 'NO')
row.setValue("Fail", 'YES')
row.setValue("Time", Start)
rows.insertRow(row)
logger.exception("Something has gone wrong!")
print "the try has failed!"

try:
    arcpy.FeatureClassToFeatureClass_conversion(failhandletest)
except:
    logger.exception("This works for the nested exception and fail has been logged in History_Table")
    rows = arcpy.InsertCursor(aTable)

for x in xrange(0, outputobjects):
    row = rows.newRow()
row.setValue("Success", 'NO')
row.setValue("Fail", 'YES')
row.setValue("Time", Start)
rows.insertRow(row)


texts: Unit of measurement from standard error prediction map?

texts: I am performing an ordinary kriging analysis and came out with the prediction standard error map. I would like to know the units of such map. Can anyone help?


texts:
I am getting this message (attached) that says:

  An ActiveX control on this page is not safe.  Your current security
  settings prohibit running unsafe controls on this page. As a result
  this page might not display as intended.

I have ArcGIS version 10, and installed it on Windows 7, 32 bit.

What should I do?


texts: How do I change the colour for selected elements?  I am selecting using select rectangle  and the selected elements are almost invisible.


texts: Distance between a line layer and several other line features (arcGIS 9.3 Arcview Licence)

texts:
I have a "Drain" line layer (with 150 lines) that I need to calculate the distance from lines to all objects within 120m. I've tried doing a spatial join and its great but I need more than just the closest features. I need all features within the 120 buffer. I believe that the NEAR tool would do this for me but unfortunately I don't have that licence Level.

Any direction you can give me would be a great help to me.


texts:
Continuously to How to remove layers that are not in visible df extent i try to remove layers that don't appear in the data frame with this code:

import arcpy,os
from arcpy import env

env.workspace = r"C:\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
    print mxdname
    mxd = arcpy.mapping.MapDocument(os.path.join(env.workspace, mxdname))
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
        if df.extent.disjoint(lyr.getExtent()):
            arcpy.mapping.RemoveLayer(df, lyr)
            print 'removed ' ,lyr
    mxd.save()
del mxd

but i get an error:

>>>
antiquities.mxd
ecology.mxd
land_use2  a2.mxd
removed  Bad_smell
removed  polygon_1
removed  polygon_1

Traceback (most recent call last):
  File "C:/Users/yaron.KAYAMOT/Desktop/python2.py", line 14, in <module>
if df.extent.disjoint(lyr.getExtent()):
  File "C:\Program Files     (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\arcobjects.py", line 144, in     disjoint
    return     convertArcObjectToPythonObject(self._arc_object.Disjoint(*gp_fixargs([second_geo    metry])))
ValueError: None
>>>

UPDATE:
i try GISGe advice in the code:

import arcpy,os
from arcpy import env

env.workspace = r"C:\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
    print mxdname
    mxd = arcpy.mapping.MapDocument(os.path.join(env.workspace, mxdname))
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
        if lyr.getExtent() != None:
            if df.extent.disjoint(lyr.getExtent()):
                arcpy.mapping.RemoveLayer(df, lyr)

    mxd.save()
del mxd

but no layer add been removed and i got:

antiquities.mxd
ecology.mxd
ta34b4.mxd
ta35_mirkamim.mxd
>>>


texts: Calculate max polygon length & width using mean aspect

texts:
I have a shapefile of polygons, and wish to calculate the maximum length (parallel to flow) and width (perpendicular to flow) - both straight lines. I calculated the mean aspect (i.e. flow direction) using DEM pixels for each polygon. The image below shows what I'm attempting (mw = width, ml = length).

Any ideas on where to start? I've found the 'Bearing Distance to Line' tool, which could potentially help to draw polylines for each polygon but am unsure regarding this and the rest.


texts:
We use an ArcSDE 10 database. We would like to log all events to export data (ArcObject) in ArcGIS Desktop 10.

Do you have an idea? a sample (ArcObject)?


texts: I would like to move a polygon's center to another location by specifying a latitude and longitude. I am using a EPSG:4326 - WGS 84 projection.


texts:
I tried raster stacking using the below code.I tried to read pixel values at a particular coordinate and few values are coming as zero where as when I tried reading the value before stacking there is some pixel value at that particular coordinate. I don't understand why my raster stacking is causing the issue.
from rasterio import plot
import os
import earthpy.spatial as es

    def image_stack(path, dest_file):
        band_fnames = os.listdir(path)
        band_paths = [path + fname for fname in band_fnames]
        arr, arr_meta = es.stack(band_paths, dest_file)
        return len(band_fnames)


texts: Is there a way to automatically snap loose points of lines to close by vertices with Quantum GIS?

texts: I have a vector layer of lines that were drawn with the snapping option off. I need to have them snapped but, there are hundreds of lines. Is there any solution to do this automatically or semi-automatically?


texts: Where can I download a map of the world for kids

texts:
I would like to download and plot some kind of fun (large) illustrated map of the world or Canada (with animals/ dinosaurs / points of interest, etc.) for my son's preschool cartography class. Anyone has a good link or map they can share?

Thanks


texts: Collecting Moran's I output values in ArcGIS Pro?

texts:
I am running an iterator in ModelBuilder that ultimately takes the output from each iteration (a point feature class) and runs the spatial auto-correlation tool (global Moran's I).  I have about 450 output feature classes, and what I would like to have is the feature class name, index value, z-score, and p-value from each iteration (of Moran's I) added as a new record in a table.  These values are reported in the .HTML report and the model run summary window, but I can't figure out how to extract them and add them to a table.

If I were doing this with a statistical model in R, I could simply call the model output values (index, p-value, etc.).  But Ideally I would like to do this in ModelBuilder.  I'm not opposed to writing a Python script, but before I go down that road I'm wondering if there's a more direct way to do it with existing geoprocessing or ModelBuilder tools.

As an additional note, I would also like to get the default calculated bandwith from the Kernel density tool output to a table similar to the Moran's I values.

I am running ArcGIS Pro 2.3.


texts: Clipping Lines Using Expression in QGIS?

texts:
I have a set of roads that I wish to clip. How ever I only need a segment of each road based on the state its located in. I can do this independently but seems like there should be a better way.

Ill try to explain it as best as I can:

I have road A, B AND C. I have states 1, 2 and 3.

All 3 roads run by all each of the states. However, I need to export only the part of road A within state 2, the part of road B within state 1 and the part of road C within state 3.

Hope its not too confusing. In reality i'm dealing with municipalities, not states but I feel this is a simpler way of explaining it.

Using QGIS 3.4 MADERIA


texts: How to work with shortest_path_shooting_star() in pgRouting?

texts:
It's been a while we are using PostGIS to store our Geometry information in our Postgres database. We have developed vertex-to-vertex routing ourselves using iGraph library. Yesterday I read about pgRouting and I am porting our current application to pgRouting. However, I face some issues in using shooting start path finder.

We stored edges and vertices information in tables and for finding a graph route with turn restriction, we are calling following SQL script:

SELECT  shortest_path_shooting_star(
    'SELECT e.id,
                e.first as source,
                e.second as target,
                ST_Length(e.geom) as cost,
                st_x(st_pointn(e.geom, 1)) AS x1, st_y(st_pointn(e.geom, 1)) AS y1,
                st_x(st_pointn(e.geom, 2)) AS x2, st_y(st_pointn(e.geom, 2)) AS y2,
                (not connected)::boolean::int * 1000000 AS to_cost,
                c.first AS rule
        FROM edges e LEFT OUTER JOIN edge_connections c ON e.id = c.second
        WHERE e.group_id = 0
        ORDER BY e.id',
    1209560, 1209653, true, false);

My inner subqery woud have following results:

  id    | source  | target  |   cost |  x1   |  y1   |  x2   |  y2   | to_cost |  rule
--------+---------+---------+----+------------------+------------------+---------+---------
1087134 |  926686 |  926687| 2.3299 |51.4675|35.728 |51.4675|35.7280|         |
1209706 | 1039731 | 1039870 | 4.005 |51.4082|35.7239|51.4082|35.7239| 1000000 | 1209564
...

In the results, when a constraint exists, rule and to_cost would not be NULL. When a turn restriction exists, the results would be exactly like what specified in the docs.

However, when I call the shortest_path_shooting_star function, results would not differ if a rule-and-to_cost exists or not. I have changed directed and has_reverse_cost parameters but nothing happened.

I was wondering if you guys could help me in this matter. Any advice or even a clue would be very much appreciated.


texts: Points not intersecting when they should

texts:
I'm new to postgis and I'm trying to show that a point intersects a polygon.

However, no matter how hard I try I can't.

Here is the point: 0101000020E6100000000000009E2FC4BF629EA494D7BB4940

Here is the polygon: https://mapit.mysociety.org/area/65731.wkt

I'm trying:

SELECT st_contains('0101000020E6100000000000009E2FC4BF629EA494D7BB4940',geom), name FROM table;

Where geom is the column with all the polygons.

As a note, the polygons are multipolygons in the database.

What am I doing wrong?


texts:
Following on from a previous question of mine, where I was given coords in degree-mins, I have successfully converted them to decimal degrees, thanks for the help.

I have created a CSV file from the stated coordinates, and then created a shapefile. The CRS was EPSG:4326 (WGS84).

The remainder of my mapping is displayed in TM65/Irish Grid (EPSG:29902). With the on the fly projections activated, everything is displayed properly. I need to provide this to a 3rd party.

How do I use QGIS to add a geometry colunm using the Current CRS? Using either QGIS 1.8 or 2.0.


texts: How to automate FGDC metadata generation in ArcGIS?

texts: Does anyone know of a tool or script to automate FGDC standard metadata for ArcGIS?  I would like to be able to select multiple feature classes and fill company info., dates, etc.


texts: How to perform SQL queries and get results from QGIS python console?

texts:
I'm writing a python script using QGIS API. I have to get results from postgres database table. The table has geometry field and other non-geometry fields.

I want to perform SQL queries on non-geometry fields and get results in my code. How to do this? Is there any class to do so?

I know about QgsDataSourceURI class but as far as I know it only gets results from geometry fields.


texts: Adding Conditional text between Dynamic "Value" text

texts:
I'm trying to show dynamic text in a text box for a map series in ArcGIS Pro. Specifically, when a map frame contains multiple boundaries (such as cities or counties). I would like to place an and between each attribute on the condition that another is present but I don't want anything to show if only one attribute is present. Delimiter is not an option because I will be using separate features.

Current dynamic text:
Project Name Unincorporated {<dyn type="table" property="value" mapFrame="Map Frame" mapMemberUri="CIMPATH=map_site_location_figures/okeechobee_county.xml" sql="Name = &apos;Okeechobee&apos;" field="Name" delimiter="and"/><dyn type="table" property="value" mapFrame="Map Frame" mapMemberUri="CIMPATH=map_site_location_figures/st_lucie_county.xml" sql="NAME = &apos;St. Lucie&apos;" field="NAME" delimiter="and"/>~ County, Florida}


texts: Why are only some tables shown in QGIS when I connect to database?

texts:
For a project I started to work with pgAdmin and QGIS, so I'm quite new in these things and I have a problem with the connection of the database to the QGIS.

I can connect myself to the database and can see the SCHEMAS "geo" and "public" and also some of the tables .. but only some of them, the one I want to connect myself to, doesn't pop up and I don't know why..
I created it like this:

-- Table: geo.casa

-- DROP TABLE geo.casa;

CREATE TABLE geo.casa
(
  id_casa serial NOT NULL,
  id_comunidad integer,
  geom geo.geometry(Point,4326),
  CONSTRAINT casa_pk PRIMARY KEY (id_casa),
  CONSTRAINT casa_comunidad_fk FOREIGN KEY (id_comunidad)
      REFERENCES geo.comunidad (id_comunidad) MATCH SIMPLE
      ON UPDATE NO ACTION ON DELETE NO ACTION
)
WITH (
  OIDS=FALSE
);
ALTER TABLE geo.casaa
  OWNER TO postgres;

I'm not sure what kind of information you need, so just tell me if you need anything else!


texts: Unable to do Map.addLayer() for all images in an ImageCollection in GEE

texts:
I found this code
Add/display all images of mycollection in google earth engine
and tried it, but the output displays four layers each with one band called constant.

Also, can you please tell me what is to be given as image.id - do I have to set any property in it?

My code is below. Note that collectYear is my image collection. How should I proceed?

geometry = /* color: #d63000 */ee.Geometry.Polygon(
    [[[76.15934177863937, 9.45738504612055],
      [76.85697361457687, 9.45738504612055],
      [76.85148045051437, 9.993390779673321],
      [76.14011570442062, 9.990685858146769]]]);
function maskS2clouds(image) {
var qa = image.select('QA60');

// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;

// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
  .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

return image.updateMask(mask).divide(10000);
}

var dataset = ee.ImageCollection('COPERNICUS/S2');
var years = ee.List.sequence(2016, 2018);
print (years);

var collectYear = ee.ImageCollection(years
.map(function(y) {
var start = ee.Date.fromYMD(y, 1, 1);
var end = start.advance(12, 'month');
var d2 =  dataset.filterDate(start, end).filterBounds(geometry)
              // Pre-filter to get less cloudy granules.
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
              .map(maskS2clouds).first();

//Applying clip to image
var ls_clipped = ee.Image(d2).clip(geometry);
var resamp_img = ls_clipped.resample('bilinear');
var B2_B11 = resamp_img.select('B2').divide(resamp_img.select('B11'));
var ratio_image = resamp_img.addBands(B2_B11);
return ratio_image;
}));
print(collectYear);

function addImage(image) { // display each image in collection
var id = image.id;
var image = ee.Image(image.id)
Map.addLayer(image);
}

collectYear.evaluate(function(collectYear) {  // use map on client-side
collectYear.features.map(addImage);
});


texts: QGIS 3.4.11: Viewing issues: Categorised symbology in point layer expands Legend in "Manage Layers Toolbar"

texts:
When I apply categorised symbology to the point layer the legend in Manager Layers Toolbar view expands to the point I cannot see any other layers in the project. This doesn't happen with other types of styling. What could be the cause of this? Is there a way to prevent this behaviour?

Edit: I use QGIS 3.4.11 version - this issue didn't happen in the version 2.18 LTR used previously
I use Simple marker point style with size 1.6 mm. Please see the picture attached for symbology settings.


texts: Enabling snapping on features created with PyQGIS

texts:
I have implemented a custom digitizing tool derivated from QgsMapToolDigitizeFeature class.
Here is the code written to save the feature into the layer:
layer.dataProvider().addFeatures([feature])
layer.commitChanges()

layer.triggerRepaint()
self.canvas().refresh()

I use the standard QGIS tool to create a new feature A. The feature A is properly added to the layer.
I use my custom tool to create a new feature B. The feature B is properly added to the layer.
Then, I use the QGIS standard tool again. Snapping works on A vertices but not on B.
If I save the project, close it and open it again, snapping works as expected on all the features.

How to make the snapping work on newly created features?
Any ideas?


texts:
I am trying to get a selection of points in my dataset, by selecting a polygon and returning all points inside it. I am using an OGC fiter, and currently I am testing it directly in the OpenLayers view of the data on geoserver.

When I use the construct:

<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:Intersects>
        <ogc:PropertyName>Geom</ogc:PropertyName>
        <gml:Envelope xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:900913">
            <gml:lowerCorner>-161672.38316912 6970834.4247659</gml:lowerCorner>
            <gml:upperCorner>-161290.19802775 6971216.6099073</gml:upperCorner>
        </gml:Envelope>
    </ogc:Intersects>
</ogc:Filter>

only the points inside the box are returned, as expected. However, I want to use a polygon, so as a test I constructed the same box with the gml:

<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:Intersects>
        <ogc:PropertyName>Geom</ogc:PropertyName>
        <gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:900913">
            <gml:outerBoundaryIs>
                <gml:LinearRing>
                    <gml:coordinates decimal="." cs="," ts=" ">
                        -161672.38316912,6970834.4247659 -161290.19802775,6970834.4247659 -161290.19802775,6971216.6099073 -161672.38316912,,6971216.6099073 -161672.38316912,6970834.4247659
                    </gml:coordinates>
                </gml:LinearRing>
            </gml:outerBoundaryIs>
        </gml:Polygon>
    </ogc:Intersects>
</ogc:Filter>

When I use this filter, however, no points are returned (neither points inside the polygon nor outside). The eventual polygon I want to use is more complex than this but I thought I'd start with a simple one. Does anyone know why this is not returning any points?


texts:
Does anyone knows a library/technique in R where you can plot a trajectory using multiple timesteps (points) where each point does have a lat/lon coordinate?
I have this data

x <- read.table("https://forms.naturwissenschaften.ch/imilast/_ERAinterim_1.5_1979_MTEX/ERAinterim_1.5_NH_M02_19790101_20121231_MTEX.txt?_ga=2.18919096.1825595846.1546710263-1112023567.1546710263", sep="", fill = T, nrows = 1000,
                     header = F, skip = 2) # download only the first 1000 rows
colnames(x) <- c("Code","CycloneNo","StepNo","DateI10","Year","Month","Day","Time","LongE","LatN","Intensity1","Intensity2","Intensity3")

Until now I was using a library called esd (developed by Meteo Norway) - which is very nice (github/metno/esd), where there was a function (read.imilast and plot.trajectory), however, that is not the greatest visualization tool). Also, I tried using openair (trajPlot) but that was quite hard to optimize for my data.

I want to plot each trajectory (column CycloneNo) under a different color if possible and more important, to have like a start and end point for each trajectory.

Any ideas?


texts: Startup() confused while converting class:IExtension into class:BaseToolbar

texts:
I developed an extension with ArcGIS 10.4.1 that I can use through the extension menu and now I need it as a toolbar. The problem is that ArcGIS launches only the startup() event from the extension and not my newly added BaseToolBar.startup()

/// <summary>
/// Classe d'extension pour le contexte
/// </summary>
[Guid("4600d1fb-cac2-4f6d-a726-261abe40dfad")]
[ClassInterface(ClassInterfaceType.None)]
[ProgId("Mrn.Taktik.ControleGis.ContexteExtension")]
public class ContexteExtension : IExtension, IExtensionConfig, IPersistVariant
{
       .  .  .  Parameters and Access variables  .  .  .

    public void Startup(ref object initializationData)
    {
        this.application = initializationData as IApplication;
        this.appStatus = this.application as IApplicationStatus;

        if (this.application != null)
        {
            contexteTakTik = ContexteTakTik.Instance;
            contexteTakTik.Initialiser(new IntPtr(this.application.hWnd));
            //// Wireup the events
            this.SetupEvents();
        }
    }

       .  .  .  #region COM Registration Functions  .  .  .
}

At this point, I simply created a BaseToolbar and copied what was inside my class into the new one. The only way to to amke it work right now is if I keep included the extension class as it is the only way for my code to enter Startup(). I tried to look online where it is first called but couldnt find anything. I need my code to go through BarreOutils.Startup() and not ContexteExtension.Startup().

/// <summary>
/// Description for ArcGISToolbar1.
/// </summary>
[Guid("bfeb9989-0c45-4db4-b12d-310efaa5d787")]
[ClassInterface(ClassInterfaceType.None)]
[ProgId("Mrn.Taktik.ControleGis.ArcGISToolbar1")]
public sealed class BarreOutils : BaseToolbar
{
.  .  .  Parameters and Access variables  .  .  .

    public void Startup(ref object initializationData)
    {
        this.application = initializationData as IApplication;
        this.appStatus = this.application as IApplicationStatus;

        if (this.application != null)
        {
            contexteTakTik = ContexteTakTik.Instance;
            contexteTakTik.Initialiser(new IntPtr(this.application.hWnd));
            //// Wireup the events
            this.SetupEvents();
        }
    }

       .  .  .  #region COM Registration Functions  .  .  .
}


texts: Add a polygon to PostGis from Google Maps

texts:
i have some experience with Google maps api, I have an app and currently calculate nearest points to locations via distance functions (great circle) and having a lot of waypoints on MySQL.

I want to start with PostGis in order to do this a better way or have been able to give a user the ability to build a polygon in GMAPS (point and click on map, saving each lat,lon, connecting the last with the first.

I've just downloaded postgresplus to start with this but don't have reference to this:
1. create a table of polygons
2. insert via php each polygon to PostGis
3. be able to query the polygon to see if a point is inside or outside a polygon

I know all this is very efficent with PostGis but have 0 experience with it.

does someone knows a link or site to learn all this basics of postgis

thanks

NA


texts:
When browsing the web and documentation I see the abbreviation lon and lng used for the longitude part of geometry.

I was wondering if there is any difference between lng and lon.


texts:
I have to work with the Shapefiles, and all is fine except data in *.dbf

When I have first deserialized file, I have thought, that I'm doing smth wrong, but when any soft like (QGIS, SAGA, ESRI shapefiles reader, GRASS GIS are reading in the same result, I have thought that source file was in the incorrect format from the beggining).

Any soft shows the next encoding from *.dbf data:

http://s15.postimg.org/lji3hjywr/image.png

Then, I've read about dBase, which is used in Shapefiles and the default encoding (LATIN1), but the docs are very confusing, because somewhere I've read about LATIN1, but at other website there is info, that for Russian maps the default encoding is: OEM866.

I've tried to reencode, but again no success.

I've downloaded the shapefiles from different sources (from bbbike and other web repositories of maps), but the same result from different resources.

My aim is, to work normally with the Moscow, Russia map, which is available here:
http://download.bbbike.org/osm/bbbike/Moscow/


texts: Is there a way to prevent losing metadata when generating a spatial view in ArcSDE 10.1?

texts:
When generating a spatial view in ArcSDE using a script, any previous metadata is dropped when the old view is dropped. Thus, any edits to metadata on the view, are lost.

Is there a way to keep the metadata around, or perhaps add the importing of metadata to view script?

I'm using Oracle as DBMS. ArcSDE 10.1 and ArcGIS for Desktop 10.2.


texts: Reading raster in Python to Clip USGS 3DEM?

texts:
I want to clip the elevation data (3DEM) downloaded from USGS website. The raster-clip tool "Clip (data management)" in ArcGIS Pro works very well. However, I cannot clip it using the python function "arcpy.Clip_management()" .

The error report that the raster file is not exist or not support. I have tried different ways to read it, such as:

~grdn39w075.tif/Band_1   This is used by the Clip tool in Arcgis Pro
~grdn39w075.tif_Band_1
~grdn39w075/Band_1

None of them can read the 3DEM file successfully.

Python code:

in_raster= "C:/Users/Desktop/Elevation/Tiles/grdn42w075.tif/Band_1"
for i in range(0,9):
 buffer = shp_list[i]
 out_id=num_buffer=re.findall('\d+', shp_list[i])
 out_path="C:/Users/Desktop/LUR_DM/GIS/Elevation/n42w75/"
 arcpy.Clip_management(in_raster,buffer,"n42w75_"+ out_id[1],
"C:/Users/Desktop/Documents/ArcGIS/Projects/LUR/LUR.gdb","99999","ClippingGeometry", "MAINTAIN_EXTENT")


texts: How to access output parameters from a Python Toolbox Tool called by .NET?

texts:
I am trying to access an output parameter from a Python toolbox tool and cannot find any information about how this is done:

I am developing an Esri Add-In for ArcGIS 10.1 SP1 for Desktop (Build 3143) and programming against ArcObjects SDK for the Microsoft .NET Framework 10.1 SP1 using C#.

I am calling a Python Toolbox tool using the following code:

// Call scripting tool from toolbox and show GUI
string assemblyFolder = Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location);
string toolbox = Path.Combine(Path.Combine(assemblyFolder, "Toolbox"), @"Amec Foster Wheeler Groundwater Tools.pyt");
string toolname = "CreateBlankGroundwaterGeodatabase";
// Open Dialog for entering parameters
IGPToolCommandHelper2 toolHelper = new GPToolCommandHelperClass() as IGPToolCommandHelper2;
toolHelper.SetToolByName(toolbox, toolname);
IGPUtilities gpUtils = new GPUtilities();
IGPTool tool = toolHelper.Tool;
IArray parameters = tool.ParameterInfo;
IGPMessages gpMessages = new GPMessagesClass();
bool pok;
toolHelper.InvokeModal(0, parameters, out pok, out gpMessages);

// What needs to be done here to get the values of the parameters, the user entered during the execution?
// How can I retrieve output parameter values here?

This is working as expected. The tool dialog opens up in ArcMap and the tool is running without any errors.

In the Python tool itself I am defining one parameter as output parameter

def getParameterInfo(self):
# [..]
param5 = arcpy.Parameter(
    displayName="out_fcpath",
    name="out_fcpath",
    datatype="GPString",
    parameterType="Derived",
    direction="Output"
)

params = [param0, param1, param2, param3, param4, param5]
return params

And in the execute function I am setting this parameter to a string value.

def execute(self, parameters, messages):
    # [..]
    arcpy.SetParameterAsText(5,"Test")

    return

But how can I retrieve that parameter by code in C# after the tool has been executed? Perhaps you can point me into the right direction?

References:

Get results from programmatically-launched geoprocessing tool
Running a geoprocessing tool using background geoprocessing with custom tool


texts:
I have a line layer with multiple traffic lines: bus, train and tram. Those lines overlaps frecuently. Every line has its own speed limit and therefore time/cost.
My goal is to do a time/cost map (isochrones) from a existing point. Previously i successfully did it with the road network for car and bike traffic, which doesn't overlap.

The tool i'm using is "v.iso.net" with 5 minutes interval, but the result differs from the previous analisys i've run, so not every line starts from the existing Point: it's not logical that a line close to the Point shows the categories equivalent to 10 or 15 minutes (blue and green), and only one has the category of 5 minutes (yellow)

To solve the problem i've tried splitting the lines every 75 meters, in case the tool was considering the whole lines, but the result didn't differed.

My question would be how to run "v.iso.net" considering all the traffic lines starts from the existing Point, and the possibility to change from one line to another in any intersecction, so it would show the minimum time/cost scenario.


texts:
I have a polyline shapefile and I am applying a Buffer tool in ArcMap 10.5 version with line_end_type='FLAT'.

The buffer shape obtained has not all the end lines with a FLAT end, some of them are round ending [Picture].

I think that the problem is where are the start and end point of lines.
Is there a process to obtained all extreme ends of the lines with a flat end?


texts:
I am quite new to QGIS but used some ArcGIS before, so I thought this should be very easy to do, but I dont get good results.

I did the following steps both in 3.8 and 3.4:
Import Raster-Layer and Shapefile, go to Raster -> Extraction -> Extract by Layermask.

Now I get a result, but both bands of the results have only 0 values.
Can someone help me pls?
Ty,
Felix


texts: Draw polygon on OpenLayers (WMS) and download data

texts:
The raster files are loaded into a PostGIS Database and connected to a WMS server, is there a good methodology to approach wherein a user can draw an AOI (Area of Interest i.e. a polygon) and download its attribute/.jpg/.img etc..

Basically similar to the USGS EarthExplorer, a small flowchart to easy the understanding, please correct me if I'am wrong in any technical details.

Raster--> PostGIS<--> WMS Server <--> OpenLayers/Leaflet <--> AOI --> Data Download


texts: QGIS 3.4.1 Extensions Repository not available

texts:
In my QGIS 3.4.1 installation (Windows 10) the extensions repository says "unavailable", like in this older post: Why I can't install plugins?.

I already added the repository with "HTTP" instead of "HTTPS", but it shows the same error.

I changed nothing on my network installation, so a proxy-problem should not be the reason.

In my older QGIS 2.18.16 installation, the repository works just fine.

Maybe anyone has a suggestion how to fix the problem?


texts: GDAL "Simple Tiling" (Or yet another - how do I tile question)

texts:
Before I start, if anyone considers this a duplicate, then that's fair enough, I'll not object, I do happen to think however that it's a question that I don't see a clear answer for either here or elsewhere.

What I need to do is to create a simple, flat set of tiles from a single large raster.

I'm aware that we have "gdal2tiles" in the os4gw toolset, but this insists on creating the simple features recommended folder structure which is not what I need to do.

I'm also aware that there are other tools (GDAL warp I think does it) that have the "Projection Window" option, allowing to extract only a given region, but for that to work I need to provide the geo co-ordinates on the command line, which means creating loops and increment things correctly, not a task that anyone working on windows really relishes :-)

I'm prepared to write my own tool if I have to, using either C and/or C# and the GDAL Libraries, but before I go down this route, I thought I'd ask here.

Put simply, is there a tool, or an option in the existing tools I may have missed that allows me to say:

c:>tilethisimage size=256x256 bigmap.tif

and end up with

tile-0-0.tif
tile-1-0.tif
tile-2-0.tif
tile-3-0.tif
tile-4-0.tif
.
.
.
.
tile-0-20.tif
tile-1-20.tif
tile-2-20.tif
tile-3-20.tif
tile-4-20.tif

Or similar, in the current directory, without the complication of having to sort through folders, figure out zoom levels or anything else?

All I need is the larger image cut into tiles of the requested size, and with the correct geo-spatial reference data embedded in them.  I don't need KML's, Virtual raster's, tab files, world files or anything else.

If it wasn't for the fact I needed to preserve the geo-spatial co-ordinate info and have those co-ordinates mapped correctly for the tile sizes, I would quite simply just use imagemagik to cut the image.


texts: What kind of raster mosaicking to use in land cover analysis

texts: I have to mosaic two rasters (Landsat ETM+ scenes) to one single image and then analyse its land cover. Which option to use to avoid seeing a distinct line in the place where two different scenes/rasters are joined? I am using ArcGIS 10.0 and eCognition.


texts: Creating Field Mapping Parameter in Arcpy Script

texts:
I'm trying to create a tool that compares two data sets. I do not want the Feature Compare tool. Most of the field names in the two data sets are the same, but not all. I want a way to map one field onto another, so I thought a Field Mappings parameter would be best suited for associating fields together. My expectation was that it would populate similar to how setting the 'Obtained From' sub-parameter to a layer populates the fields list with all available fields from that layer, but that did not happen. I want functionality almost exactly like the Append tool. It populates the 'Field Map (optional)' window with a list of field maps where the user can add input fields.

I've tried working with the tool validator class, in the UpdateParameters method, I've added code to create a FieldMappings object and fill it with FieldMap objects created from the fields in another parameter, but I can't get the window on the tool itself to fill with the FieldMap objects. I've debugged to the point where I know for a fact the FieldMappings object is filled how I want it. I just can't get it to display.

Here is my code for the updateParameters method. param[0] and param[1] are feature layers. param[2] is a field mappings parameter.

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    if self.params[0].altered:
    self.params[2].value = arcpy.FieldMappings()
    for field in [a for a in arcpy.ListFields(self.params[0].value.dataSource) if not a.required]:
        fm = arcpy.FieldMap()
        fm.addInputField(self.params[0].value.dataSource,field.name)
        self.params[2].value.addFieldMap(fm)
    return

I have also tried this:

def updateParameters(self):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        if self.params[0].altered:
            fms = arcpy.FieldMappings()
            for field in [a for a in arcpy.ListFields(self.params[0].value.dataSource) if not a.required]:
                fm = arcpy.FieldMap()
                fm.addInputField(self.params[0].value.dataSource,field.name)
                fms.addFieldMap(fm)
            self.params[2].value = fms
        return

I can't find a single post on creating a field mapping parameter in a script similar to the append or spatial join tool.

If the autopopulation is not possible, even getting the FieldMapping parameter to work properly would be an acceptable middle ground. Right now, the 'Add Input Field' never appears when right clicking a field map parameter. If I set the Field Mapping parameter to single value, clicking the window gives 2 options, Add Output Field and Reset. If I add an output field, right clicking that gives more options, Delete, Rename, Merge Rule, and Properties. None of which allow me to add an input field to the field map.


texts: There is no overlap with watershed grid in HRU analysis?

texts:
I have to calibrate a SWAT model. I have my watershed, land-use and soil data file in projected coordinate system. The coordinate system is NAD_1983_UTM_16N. When I upload the land-cover data on HRU analysis tab, an error appears:

  There is no overlap with watershed grid.

I have looked up on other resources but could not find a solution.

Any solution?

Here is the picture of my watershed and land-use layers which are perfectly overlapping:


texts: Leaflet's map.project(latlng, zoom) without a leaflet map instance

texts:
I am working with a leaflet app and I'm shifting a lot of the GIS calculations to the back end (a nodejs app). I need to be able to get the RGB value of a pixel given the latlng and the map zoom.  (I'm using the Mapbox RGB Encoded Dem and I need elevation on a pixel-by-pixel basis).  In order to do that, I need the XY value of the pixel relative to the map origin.  This is conveniently done using leaflet's map.project function.   You can see a quick example of this here: codesandbox example.
Now I want to be able to get these XY coordinates on the back end, and I won't have the leaflet map instance to help me.  (I could pass the L.map instance to the back end, but I feel this is overkill.)
How can I write a function from scratch (or with some small / simple GIS libraries) to calculate the XY pixel value from a given latlng and zoom?  I feel this is probably a simple or even duplicate question but I'm not finding any coord transformation discussions that account for map zoom, or that are specific to a latlng's associated pixel on a webmap for a given projection and xyz coord.


texts:
I try to fit my shape to the google map.

require(rgdal)

in.dir <- getwd()

    sh <- readOGR(in.dir, layer = "Ytor",
p4s = "+proj=tmerc +lat_0=0 +lon_0=15.80827777777778 +k=1 +x_0=1500000 +y_0=0 +ellps=bessel +units=m +no_defs")

sh@proj4string
sh2 <- spTransform(sh, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs"))

sh2@proj4string

 sh2.df <- fortify(sh2)

library(ggplot2)
library(ggmap)

swe <- get_map(location = c(15.3, 63.63, 15.5, 63.68),  maptype = 'satellite')
ggmap(swe) + geom_polygon(aes(x = as.numeric(long), y = as.numeric(lat), group = group), data = df, color = "white", fill = "white", alpha = 0.1)

Finally I got the map as below

The borders of the shape are not going as the ggmap borders. What did I wrong that the borders does not fit?


texts: Is it possible to export styled vector data as tif in QGIS without losing transparency?

texts:
I have vector point data which I have styled using svg images, see below with the globe in the background:

And without:

I want to export the image without the globe as a tif with transparency. The white should be 100% transparent while the edges around the data should keep there semi-transparency.

I have tried "save image as", which doesn't keep transparency and is to low resolution. Then I tried the following command:

.\qgis.bat --project D:\test\hurricanes\hurricanes_test.qgs --snapshot D:\test\hurricanes\hurricanes_new.tif --width 30000 --height 15000

I was able to increase the resolution this way, but again it doesn't keep transparency. I have tried using the background setting in the print composer, but it doesn't keep the semi-transparent edges around the data. Does anyone have any ideas?


texts: Showing QgsComposition created in QGIS Python console

texts:
I'm trying to develop a similar script to How to make serial maps from template?. I'd like to modify elements in the composition programmatically for each map that I want to print. In developing this I'd like to see the elements I'm modifying in the console without printing off the map. I was expecting to see a new print composer window open when I ran the below code from the QGIS console in 2.18.2, but nothing happens (no error message appears either).
from qgis.PyQt.QtXml import QDomDocument
def loadPrintComposerTemplate(template):
    '''Load a print composer template from provided filename argument

    Args:
        template: readable .qpt template filename

    Returns:
        myComposition: a QgsComposition loaded from the provided template
        mapSettings: a QgsMapSettings object associated with myComposition'''
    mapSettings = QgsMapSettings()
    myComposition = QgsComposition(mapSettings)
    # Load template from filename
    with open(template, 'r') as templateFile:
        myTemplateContent = templateFile.read()

    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent)
    return myComposition.loadFromTemplate(myDocument), mapSettings
composer, mapSet = loadPrintComposerTemplate(template)

Is there a way to make the print composer appear?


texts:
How can set zoom level scale for tile map?

I want set zoom level scale like google map.

Set zoom level according to tiles.


texts: When I save the shapefiles and the project the shapefiles do not save so I loose a lot of work every day. The problem starts when qgis starts to delay the save of the layers after long time in edit mode, so at the end of the day I save the project and the layers and the next day I don't have the last information in the layers that I saved. I have also set up the auto save for the 30 minutes but still doesn't save. What is the problem? I think it is the cache information. I am using qgis is 2.18.18.


texts: PostGIS via GeoDjango: function ST_DistanceSphere(geometry, geometry) does not exist

texts:
I've been using PostGIS via GeoDjango (Django v3.0.7), happily until just now.  I have a query (Django .filter()) that looks like this:
.filter(point__distance_lte=(dlPt, D(m=CloseRadius)))

when it invokes the PostgreSQL query I get this error:
django.db.utils.ProgrammingError: function ST_DistanceSphere(geometry, geometry) does not exist

dlPt is indeed a Point, <class 'django.contrib.gis.geos.point.Point'>, and so are the point attributes in the objects being searched.
Any ideas what the issue might be?


texts: Representation of Summary statistics in QGIS composer

texts:
Is it possible to display a summary of statistics from a certain attribute within a feature class within the print composer using an sql statement. I have an shapefile "projects" with a number of features and attributes such as location, author, value, etc etc. I would like to provide a simple table or text box within the composer that displays statistics based on a certain attribute. For example i would like to show the sum of all the "value" attributes for a given "author". something like

Sum('projects'.'value') where 'author'=XXX

Therefore, if i change some of the attributes in the table, such as the value of a project or author of a project, the text in the composer would update automatically.

Ideally i would like a composer that shows on the map all of the features (coloured by author for example) and then a summary table for each author with a feature count, total value, and any other summaries etc.


texts: Getting overlap of supply catchments using ArcMap

texts:
I need to know, how much each of my calculated catchments is overlapped by the others, to assess the supply of the commodity per source point.
I calculated supply catchments for a total of 110 point features derived from an accessibility analysis for a specific commodity. That makes 110 catchment rasters with 3 classes OR 330 feature classes if converted to polygon by using the 3 classes.
I know the methodology to do it with intersect, but with this amount of layers it just takes too much time.
Any idea how to do it?


texts: Cutting raster using ArcGIS for Desktop?

texts:
I would like to cut a raster file that I have.

Here I attach a example of how it is original:

And this is how I would like it would be. (I have cut the image on Paint):

Can someone give me a hint on how can it be done? Is it simple like a "select" and delete?


texts: Nested if statement in QGIS label expression

texts:
I'm trying to create labels for features, based on the value (or missing value) of three different columns and combinations can be possible.
I have three columns: Veg_dom, Veg_codom and Veg_onderg.
When all three columns have a value for a feature, I want the label to be as follows:
Veg_dom$value + Veg_codom$value + (Veg_onderg$value).
Other possibilities should be
Veg_dom$value + Veg_codom$value & Veg_dom$value + (Veg_onderg$value).
I use the following expression in the label expression window:

The labels however show as follows:

The only one correct is the first statement if all three columns have values.
Am I doing something wrong? Are nested IF statements allowed?
Expression:
if( Veg_codom <> ' ' AND  Veg_onderg <> ' ',
Veg_dom + '+' + Veg_codom + '+' + concat('(', Veg_onderg, ')'),

if( Veg_codom = ' ' AND  Veg_onderg <> ' ',
Veg_dom + '+' +concat('(', Veg_onderg, ')'),

if( Veg_codom <> ' ' AND  Veg_onderg = ' ',
Veg_dom + '+' + Veg_codom,

Veg_dom)))

Cels hold character-data:


texts: How can I calculate averaged response curves for Maxent models fitted using maxnet package in R?

texts:
Could someone please point me towards an efficient approach?
plot and response.plot (maxnet package) can produce response curves for individual models, but feeding a list of model replicates, as generated when fitting models in a for-loops using different sets of occurrences, fails.
Any idea or comments?

My models are created in a similar way as in the following example:

require(maxnet)
data("bradypus")
presences <- which(x = bradypus$presence == 1)
absences <- which(x = bradypus$presence == 0)

my_models <- list()
set.seed(123)
for (i in 1:3) {

maxnet_data <- bradypus[c(sample(x = presences,
                                 size = 25,
                                 replace = F),
                          sample(x = absences,
                                 size = 25,
                                 replace = F)),]

my_models[[i]] <- maxnet(p = maxnet_data[, 1], data = maxnet_data[, -1],
                         f = maxnet.formula(p = maxnet_data[, 1],
                         data = maxnet_data[, -1], classes = "lq"))
}


texts: Creating relationship class in QGIS for uploading to ArcGIS Online?

texts:
Since I don't have the advanced license for ArcGIS Pro, I'm trying to create relationship class in other programs.

Would ArcGIS Online support relationship (O-M) created through another program (such as QGIS or FME)?


texts: window.open and adding multiple maps to DOM of new window

texts:
I have created an HTML template to open from my existing arcgis JavaScript application. The HTML template is to be populated with several maps and text based on a user selection of the main app. My question is, what is the best way to go about copying the same graphics and basemaps from the main application to the new window? My current code adds a new map and I am able to add the basemap as a layer, but it is removed in the main application. I also have not been able to add the graphic layer, although i get no errors when I add to the maps of the new windows.

here is my current code:

  var center_left = (screen.width / 2) - (1100 / 2);
  var center_top = (screen.height / 2) - (600 / 2);
  var adminWindow = window.open('', '_blank', 'width=1100,height=600,menubar=yes,scrollbars=yes,left=' + center_left + ',top=' + center_top);
  var adminDoc = adminWindow.document;
  var adminBody = adminDoc.body;
  adminDoc.write(adminReportTemplate);
  adminDoc.close();

  var parcelTable = this.getParcelTable();
  adminDoc.getElementById('parcelTable').innerHTML = parcelTable;

  var parcelExtent = this.results.features[0].geometry.getExtent();

  var basemapLayer = this.map.getLayer(this.map.layerIds[0]);

  var parcelGraphic = this.map.getLayer(this.map.graphicsLayerIds[8]);

  var parcelMap = new Map(adminDoc.getElementById('parcelMap'), {
    extent: parcelExtent,
    basemap:basemapLayer  // Error: Map.setBasemap: Unable to add basemap: "".
  });
  parcelMap.addLayer(basemapLayer); //works but removes basemap from this.map
  parcelMap.graphics.add(parcelGraphic); //no graphics visible, no error


texts: How to apply a rectangular buffer on spatial polygons in R/qGIS

texts:
I am trying to find a way to apply a rectangular buffer on my spatial polygon dataframe in R. At the end, I would like to be able to shrink my polygons by specified units on both sides. So far, I have tried gBuffer function from rgeos package. However the issue with gbuffer is that it only takes single argument of width and hence applies a uniform buffer on all sides (like buffer by percentage in qgis). My original polygon file has 1800 shapes. Here is a small dummy set to illustrate the problem in R:

library(rgeos)
# simple example, from scratch:
Sr1 = Polygon(cbind(c(1,4,4,1,1),c(1,1,4,4,1)))
Sr2 = Polygon(cbind(c(1,4,4,1,1),c(6,6,9,9,6)))
Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
SpP = SpatialPolygons(list(Srs1,Srs2), 1:2)
#left panel on the figure
plot(SpP)
plot(gBuffer(SpP, byid=TRUE, width=c(-0.3, -0.3)), add=T)

Here is the output figure:

PS: I am also open to qgis solutions.


texts: mapCanvas.saveAsImage() doesn't work properly within loop in standalone script (only .PGW)

texts:
I have a problem with the export of QGIS 3 interface screenshots. I run stand-alone scripts and one of them launches QGIS and runs script.py with the following code :

mapCanvas = iface.mapCanvas()

selectedfeatures= layer.getFeatures(QgsFeatureRequest(request))

for feat in selectedfeatures:
    layer.selectByIds([feat.id()])
    mapCanvas.zoomToSelected(layer)
    mapCanvas.refresh()
    mapCanvas.saveAsImage(path)

with the command line :
qgis-bin.exe --nologo --project path_project --code script.py

But the images are not always saved (random files are saved, sometimes nothing), only the .pgw like if the interface had no time to refresh... I already read solutions like with refresh SLOTS, or QTimer but it doesn't work for me.


texts: What is a general purpose of "Update SQL layer" in QGIS 3?

texts:
Can someone explain to me, what is a general purpose of "Update SQL layer"?

When I "update" it becomes not editable... (see gif ▼)