texts: Select by expression only for selected features using PyQGIS

texts:
I'm facing a nested problem using PyQGIS (for QGIS master, but it can be replicated in any version of QGIS 2).
I have a vector layer with some selected features and with a QgsFeatureRequest I can easily loop only on the selected features:
vl = iface.activeLayer()
# get the list of selected ids
ids = vl.selectedFeatureIds()

# create the request with the selected ids
request = QgsFeatureRequest()
request.setFilterFids(ids)

# use only the selected features in getFeatures
for i in vl.getFeatures(request):
    print(i.attributes())

Then I build an expression based on a field of the attribute table and I can use QgsFeatureRequest() to filter only the features that match that expression:
exp = ''' "Piper" = 1 '''
request = QgsFeatureRequest()
request.setFilterExpression(exp)

for i in vl.getFeatures(request):
    print(i.attributes())

but what if I want to filter with the same expression only the already selected features?
So actually use both requests in the same getFeatures loop?


texts: Only integer part of field values read by pyshp when creating point shapefile from QGIS

texts:
While using flopy that bases some of the functions in the export.shapefile_utils module on pyshp (namely shp2recarray), I've come across a behavior I cannot understand.
Whenever I read some point data shapefile, the records that I can later inspect in Python only show integer values.
I've verified that when writing point shapefiles with pyshp, the Writer class writes only integers values unless the decimal field is specified. Thus, of the following three blocks of code the only one that generates a shape where I can see the decimal figures in the attribute table of the resulting file, is the last one:
w = shapefile.Writer('./test/testfile1')
w.field('X', 'N')
w.field('Y', 'N')
w.field('Z', 'N')
w.point(122.5, 37.5)
w.record(122.5, 37.5, 23.5)
w.close()

w = shapefile.Writer('./test/testfile2')
w.field('X', 'F')
w.field('Y', 'F')
w.field('Z', 'F')
w.point(122.5, 37.5)
w.record(122.5, 37.5, 23.5)
w.close()

w = shapefile.Writer('./test/testfile3')
w.field('X', 'N', decimal=3)
w.field('Y', 'N', decimal=3)
w.field('Z', 'N', decimal=3)
w.point(122.5, 37.5)
w.record(122.5, 37.5, 23.5)

This seems to be consistent with pyshp documentation that states:

Reading Records
A record in a shapefile contains the attributes for each shape in the
collection of geometries. Records are stored in the dbf file. The link
between geometry and attributes is the foundation of all geographic
information systems. This critical link is implied by the order of
shapes and corresponding records in the shp geometry file and the dbf
attribute file.
The field names of a shapefile are available as soon as you read a
shapefile. You can call the "fields" attribute of the shapefile as a
Python list. Each field is a Python list with the following
information:
Field name: the name describing the data at this column index.
Field type: the type of data at this column index. Types can be:
    "C": Characters, text.
    "N": Numbers, with or without decimals.
    "F": Floats (same as "N").
    "L": Logical, for boolean True/False values.
    "D": Dates.
    "M": Memo, has no meaning within a GIS and is part of the xbase spec instead.
Field length: the length of the data found at this column index. Older GIS software may
truncate this length to 8 or 11 characters for
"Character" fields.
Decimal length: the number of decimal places found in "Number" fields.

However when creating simple point shapefiles in QGIS with real value fields, in the way depicted by the figure, I only get the integer part of values in each field read when reading the file with pyshp.

Is there any setting I'm missing when creating the shapefile in QGIS that is preventing the decimal part of the fields being read by pyshp?
I'm reading the shapefile with the estandar Reader
sf = shapefile.Reader('./test/somefile')
sf.record(0)
# Out[29]: Record #0: [122, 37, 23]


texts: Build the edge-vertex topology from a polyline feature class in ArcGIS

texts:
I have a polyline feature class called line_1. What I want to achieve is the edge_vertex topology. Specifically,

I want to generate a new polyine feature class called line_1_new, by adding two new columns to line_1, start_v_id and end_v_id.
I want to generate a point feature class from line_1. The point feature has an attribute called v_id, corresponding to the start_v_id and end_v_id attributes.

I have used the method of Feature Vertices To Points to generate a point feature class of both ends of the lines. However, the point feature does not contain the topology explicitly. For example, if two lines intersect at a junction, this method will generate two points corresponding to the junction, even though these two points have the same coordinates.


texts: Calculating centroids of polygon using ArcPy?

texts:
I am trying to calculate the centroids of a polygon and add the x- and y-coordinate to their respective fields:

arcpy.env.workspace = r"C:/..."

input_fc = "shapefiles/precincts_v.shp"
arcpy.AddField_management(input_fc, "X", "DOUBLE")
arcpy.AddField_management(input_fc, "Y", "DOUBLE")

arcpy.CalculateField_management(input_fc, "X", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
arcpy.CalculateField_management(input_fc, "Y", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[1]", "PYTHON")

but I an invalid syntax error. Where is the mistake?

Edit: exact error message

Traceback (most recent call last):
  File "C:\Users\mquentel\Dropbox\PLSS\CA voting data\Python\CONSTRUCT_merge_voting_and_census.py", line 17, in <module>
    arcpy.CalculateField_management(input_fc, "X", "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 3360, in CalculateField
    raise e
ExecuteError: ERROR 000539: SyntaxError: invalid syntax (<expression>, line 1)
Failed to execute (CalculateField).


texts: How to apply colour ramp to user-defined range of values for contour lines?

texts:
I have contoured elevation data, across a number of shapefiles, at 10m intervals. Each tile has a different range of values (min - max).
Under Style > Categorized, I can only apply the colour ramp to discrete values in the current tile. I can add further categories to extend the range min and max to encompass all of my tiles, but I don't see how to apply the colour ramp across this whole range.

So, I don't seem to be able to create a single categorization that I can apply across all of the tiles.

(Also, it's a bit annoying that the 'Value' sort doesn't recognise that these values are from an integer field, and seems to sort them as strings, so I can't re-order the categories in numerical order... off to hub.qgis.org/issues I go)

I'm using QGIS 2.6.0

any help would be gratefully received.


texts: Whats a Python Script/Method to use to get the number of fields with data in an attribute table?

texts:
I am doing some data migration and need the number of records with actual data, sans null and blank fields, without switching back and forth with select by attribute.

Ideally I would like to run a script for a table or feature class that gives me the count of each field that has data.


texts: Selecting different field in attribute table for kriging in ArcGIS for Desktop?

texts:
I have created a shapefile with several columns of data in the attribute table but I can not select any of the columns for a kriging except the 4 required to display (lable, northing, easting, elevation)

my attribute table is something like

well.id  northing  easting  elevation  data1  data2  data3  data4

I wish to krig the data in the data columns however my options for Z value field in the kriging window are only the first four column titles...

any suggestions?


texts: Failed to create the service The submitted field 'KML_FOLDER' does not exist in the table 'db_8bog9.hsu_gk2bm.company_layers_wfl1_horizon_lines'

texts:
Receiving an error when trying to share a webmap to portal using ArcGIS Pro 2.4.3,

Failed to create the service 'Hosted/Company_Layers_WFL1.FeatureServer'. The submitted field 'KML_FOLDER' does not exist in the table 'db_8bog9.hsu_gk2bm.company_layers_wfl1_horizon_lines'.

The thing I cannot figure out is the specific table does not have a field named 'KML_FOLDER" and I looked through all ~80 of the layers in this webmap and none of them have a Field with that name. My Enterprise set up is 10.7.1 and SQL Server.

Below are the fields in this layer that is being referenced in the Error. I cannot simply remove this layer as it is rather important.


texts: How to point out a city from a map image of world with longitude and latitude?

texts:
have a map image 800px wide. How can I figure out a city.

like New York has 40.75 latitude and -73.98 longitude. I have to pin a mark over this point.

Before I have done with x_axis and y_axis which gives problem on more cities figure because i manually had to input values.
now i have database with all cities with latitude and longitude.
or is there any way that i can convert latitude and longitude to x and y axis.

Thank you

EDIT

x = image.width * (longitude + 180) / (2 * 180)
y = image.height * (latitude + 180) / (2 * 180)


texts: Trying to search for gdal executable versions anywhere like the version from Unofficial Windows Binaries for Python Extension Packages before it was change to wheel but no luck(e.g. GDAL-1.10.1.win-amd64-py2.7.exe). Are gdal executables still available? Or was this replaced with wheel?


texts: Using georeferencer plugin of QGIS (gdal)?

texts:
I need to rotate the tiff image that you can see below (basically, consider that the RED rectangle belongs to the tiff: I want to rotate it so that the RED rectangle aligns with the GREEN one).

To do so I tried to use the georeferencer plugin (gdal). But when I try to align it using transformations like 'projective' or 'linear', the georeferencer acts as if it was working, but the output file cannot be readed anywhere as if it was an empty file (even though it is not).
Could the problem come from the fact that my tiff is quite big (1.7 GB)? (Because the same method worked 2 weeks ago with a smaller tiff)

Is there actually a plugin,a script to run in the python console, or any other (free) software that could help me rotating my tiff?

As I have many tiff to rotate, it would be much faster.


texts: How to create perpendicular line from point on one line to an intersection with the nearest polygon layer?

texts:
I am trying to create laterals from a layer of sewer mainlines and have yet to Google a thoroughly helpful answer.  Hoping to find some help here.

I have a layer of mainlines that have a direction associated with them pointing from upstream to downstream.  I have tabular data that represents the intersection of a lateral line with the mainline.  It has the ID of the mainline, the distance from the upstream manhole where it intersects and the side of the mainline that it is on (numbers 0-12, directly downstream being 0/12 or directly upstream being 6, 7-11 being to the left of the downstream direction of the line).  I need to create perpendicular lines from these intersection points to the nearest polygon in my parcel layer and extend that lateral line 2 feet into the parcel.  I also need to add that parcel number to the lateral in my attribute table and label the laterals ID_1, ID_2, ID_3, etc. on each mainline.

So far in ArcMap 10.0 sp4 ArcInfo license, I have created a set of points representing those intersections on the mainline.  With Linear Referencing, I created routes out of the mainlines and then made an event layer of lateral lines using the distance from the "zero" measure of the mainline.  Now I'm stuck.  How do I make perpendicular lines on the correct side of the mainline for these 45k laterals to cross 2 feet into the nearest parcel polygon?


texts: QGIS OpenLayers OSM plugin not matching with my imported GIS information

texts:
I have imported GIS data into a new QGIS project for our local city website. The files are shapefiles and inside the prj it shows:

  PROJCS["NAD_1983_StatePlane_Florida_East_FIPS_0901_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",656166.6666666665],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-81.0],PARAMETER["Scale_Factor",0.9999411764705882],PARAMETER["Latitude_Of_Origin",24.33333333333333],UNIT["Foot_US",0.3048006096012192]]

I set the CRS to NAD83 Florida East - ft on each layer I imported and set the OpenStreetMap layer the plugin created to the same CRS. The project settings also has the same CRS however the GIS data I imported is no where near where it is supposed to be. It is on the other side of the world.

What am I missing?


texts: How to deal with "ERROR: GEOSUnaryUnion: TopologyException: Self-intersection"?

texts:
I'm using Postgres 9.4 and PostGIS 2.1. My table looks like this:

------------------------------------------------------------------
 ogc_fid      | integer                | not null default
 wkb_geometry | geometry(Polygon,4326) |
 county       | character varying      |

I am trying to create the union of several thousand polygon fields:

select ST_UNION(wkb_geometry) from mytable where county='Glos';

But I get this error:

ERROR:  GEOSUnaryUnion: TopologyException: Input geom 0 is invalid:
Self-intersection at or near point -2.3936967722194602 51.540060494135503 at -2.3936967722194602 51.540060494135503

Why could this be happening? And what can I do about it? Is there a way that I can check whether the input polygons are valid?

I don't mind if the fields overlap: I just need to merge them all down (as it were) to create the overall union.

UPDATE: I'm trying this:

SELECT ogc_fid, ST_isValidReason(wkb_geometry) FROM mytable WHERE (county='Glos') AND (NOT ST_isValid(wkb_geometry));

Results:

   ogc_fid |                      st_isvalidreason
-----------+------------------------------------------------------------
  53270254 | Ring Self-intersection[-2.42234450631667 51.5302222606785]
  20669431 | Ring Self-intersection[-2.42680440846532 51.5321194786993]
  20851050 | Self-intersection[-2.39369677221946 51.5400604941355]


texts: Calculating area for polygons in arctic areas

texts:
We are evaluating areas of ecological importance in the arctic. The asset areas are of different sizes (up to 10 square kilometres) from latitudes 64 to 80 degrees north. The areas also vary in longitude from -10 degrees to +40 degrees.

Since the areas calculation are important in the evaluation of the assets we need a conform method to calculate their areas. What would be the best approach? Which projection would be better?


texts: why preview image is blank for ImageMosaicJDBC (postgis raster) layer?

texts:
I follow http://docs.geotools.org/latest/userguide/library/coverage/pgraster.html to publish a ImageMosaicJDBC layer in geoserver.
the layer seems published ok and saved, but the preview of this layer does not give me a meaningful overview image, only a blank image with the correct lat lon extent.  Also I don’t know why I got two bands when I get feature info, the geotiff that goes into postgis raster table only has one band.  In the contrast, using Quantum GIS Lisboa I can see greyscale image after adding layer for the Postgis raster table in postgres.

below is what I have done:

step 1: Translate netCDF into GeoTiff:
$ gdal_translate.exe -of GTiff -scale -32766 32766 -32766 32766 -ot INT16 gebco_08_-12_-51_11_54.nc gebco08part2.tif

step 2: Load GeoTiff into Postgres as PostGIS raster table:
C:\temp>"C:\Program Files\PostgreSQL\9.2\bin\raster2pgsql.exe" -s 4326 -I -C -M gebco08part2.tif | "C:\Program Files\PostgreSQL\9.2\bin\psql.exe" -d postgis20 -U postgres

step 3: Create master table and tile table in Postgres:
create table gebco_master2 (NAME varchar(254) not null,
         TileTable varchar(254)not null,
         minX FLOAT8,minY FLOAT8, maxX FLOAT8, maxY FLOAT8,resX FLOAT8, resY FLOAT8,
         primary key (NAME,TileTable))

insert into gebco_master2(NAME,TileTable) values ('gebco_cover','gebco08part2')

step 4: Configure mosaic JDBC plugin in geoserver:
gebco.pgraster.xml:
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ImageMosaicJDBCConfig [
      <!ENTITY mapping PUBLIC "mapping"  "mapping.pgraster.xml.inc">
      <!ENTITY connect PUBLIC "connect"  "connect.pgraster.xml.inc">
]>
<config version="1.0">
      <coverageName name="gebco_cover"/>
      <coordsys name="EPSG:4326"/>
      <scaleop  interpolation="1"/>
      &mapping;
      &connect;
</config>

connect.pgraster.xml.inc:
<connect>
  <!-- value DBCP or JNDI -->
  <dstype value="DBCP"/>
  <!--   <jndiReferenceName value=""/>  -->
  <username value="postgres" />
  <password value="password" />
  <jdbcUrl value="jdbc:postgresql://localhost:5432/postgis20" />
  <driverClassName value="org.postgresql.Driver"/>
  <maxActive value="10"/>
  <maxIdle value="0"/>
</connect>

mapping.pgraster.xml.inc:
<spatialExtension name="pgraster"/>
 <mapping>
      <masterTable name="gebco_master2" >
              <coverageNameAttribute name="name"/>
              <maxXAttribute name="maxx"/>
              <maxYAttribute name="maxy"/>
              <minXAttribute name="minx"/>
              <minYAttribute name="miny"/>
              <resXAttribute name="resx"/>
              <resYAttribute name="resy"/>
              <tileTableNameAtribute name="tiletable" />
      </masterTable>
      <tileTable>
              <blobAttributeName name="rast" />
              <keyAttributeName name="rid" />
      </tileTable>
</mapping>

step 5: Then add ImageMosaicJDBC store and layer in geoserver:
URL Connection Parameters:
file:coverages/gebco.pgraster.xml

I believe the postgis raster table in postgres database should be correct otherwise quantum gis won't be able to view it, the question is have I configured geoserver correctly, why the preview image is blank?


texts: Shapefile missing from 'Add Vector Layer' file type in QGIS?

texts:
Ran into a bizarre problem in QGIS ( 2.18.4). I am trying to load a shapefile with extension *.shp.

When I choose Add Vector Layer -> Browse, I do not see any option to load ESRI Shapefiles.

I must be missing something very elementary.


texts: Removing "(optional)" from parameter description in Python script tool dialog?

texts: I am using ArcGIS 10.2 and create a script tool in toolbox for attribute search on a layer, all the searching criteria is not compulsory, and user can combine multiple search criteria. But in the window all parameter together with (optional), I want to know whether there is a way to remove the (optional), the picture is attached.


texts: Obtaining Orientation of Points to Another Layer?

texts: I have a vector layer with points denoting receptors. I have a second vector layer (polygon) delineating a site boundary.  Utilizing QGIS v2.14, I am looking for a method to get the orientation of the points to the polygon. I.e. in degrees what is the orientation of each point to the site boundary polygon.


texts:
Somehow my postGIS template got overwritten by one with data from one of the databases I previously created.

I've tried everything I read on goggle but could not delete or replace the template with the correct one... I am hoping that someone here would know the best way to fix this problem without having to do a complete a database software.

Thanks
Chris


texts: Geoserver &transparent=true equivalent with SLD?

texts:
I have some polygons (country borders) and using the query parameter "transparent=true" in my request, I can see what is behind the layer, where no border is.
That is great!

On the other hand, I have to style the border and I should not do it on the server side. So I send an SLD with my map request.
Styling the borders, filling the polygons with color, that works fine.

Working with requests from the SLD Cookbook like this one

 <FeatureTypeStyle>
     <Rule>
       <PolygonSymbolizer>
         <Fill>
           <CssParameter name="fill">#000080</CssParameter>
           <CssParameter name="fill-opacity">0</CssParameter>
         </Fill>
         <Stroke>
           <CssParameter name="stroke">#FFFFFF</CssParameter>
           <CssParameter name="stroke-width">2</CssParameter>
         </Stroke>
       </PolygonSymbolizer>
     </Rule>
   </FeatureTypeStyle>

do not work for me. I have still the white background.
I also tried to leave  <CssParameter name="fill">#000080</CssParameter> and also the whole fill-part out. But no success.

I also tried to combine the &transparent=true parameter and the SLD. But that was also without success.

Last but not least I tried to find some configuration on the server side.
Thought maybe I can set the transparency per default, so If I want to fill the polygon I can explicitly say this with SLD.
But as you might already guess...no success.


texts:
I am trying to install the opengeo suite (community) for a project I have. The thing is that my Linux distribution is Arch Linux, for which the package manager is pacman and not apt-get as it is the case in Ubuntu.

Is it possible to install it on Arch Linux? The installation manual says that it provides packages for both debian and red hat, but is possible to build a package from sources somewhere? I know they have a git repository, however, I cannot see how to build a package from it (my knowledge of writing PKGBUILDs is slim at best...)

Thank you


texts: ArcObjects edit/update ODBC Oracle table

texts:
Can an ODBC Oracle table attached by OLE within ArcMap be edited/updated by using ArcObjects? If so, what are the procedural steps for doing this?

Are there code examples out there showing how to do this?


texts: Changing the color of specific area of shapefile using Python

texts:
I have this map and I need to emphasize with a different color some cities from this map:

This is my code:
import json
import time
import geopandas as gpd
import adjustText as aT
import matplotlib.pyplot as plt
from descartes import PolygonPatch
from shapely.geometry import LineString

borough = gpd.read_file('London_Borough_Excluding_MHW.shp')
borough = borough.to_crs({'init': 'epsg:4326'})

borough["center"] = borough["geometry"].centroid
borough_points = borough.copy()
borough_points.set_geometry("center", inplace = True)

texts = []

ax = borough.plot(column=None, cmap='ocean', edgecolor = "black", linewidth = 1, alpha=0.5,
                  categorical=False, legend=False, axes=None, figsize=(40,20))

for x, y, label in zip(borough_points.geometry.x, borough_points.geometry.y, borough_points["NAME"]):
    texts.append(plt.text(x, y, label, fontsize = 12))

aT.adjust_text(texts, force_points=0.3, force_text=0.8, expand_points=(1,1), expand_text=(1,1),
               arrowprops=dict(arrowstyle="-", color='grey', lw=0.5))

x, y = df2['lon'].values, df2['lat'].values
ax.scatter(x,y, marker="*", color='r', zorder=5, alpha=0.7, s=18)
ax.set_ylim([51.65, 51.40])
ax.set_xlim([-0.4, 0.2])

plt.show

Could someone show me the right direction to do it?


texts: QGIS: prefixPath() gives NULL in PyCharm

texts:
I have installed QGIS 2.18, QGIS 3.8 and OSGeo on a Windows machine.

I tried to set up a PyCharm environment for PyQGIS coding, but I can't initialize the OGR provider to import shapefiles.

Here's the code:

qgs = QgsApplication([], False)
qgs.setPrefixPath(r'C:/PROGRA~1/QGIS3~1.8/apps/qgis', True)
qgs.initQgis()

When I try to run

qgs = QgsApplication.setPrefixPath(r'C:/PROGRA~1/QGIS3~1.8/apps/qgis', True)
print(qgs)

It returns None. If I have installed three versions of QGIS in Windows there are some problems with prefixPath()?


texts: Compute a distance matrix in R with geosphere

texts:
I am trying to compute the a distance matrix between two sets of points in two shapefiles (S1 and S2) in R. I am new to geographical analysis in R, having previously only used QGIS and ArcMap.

I am using to do so the packages geosphere and rgdal. My two shapefiles are in the projection  EPSG:102023. When I try to compute the distances using distm, I get the error:

 Error in .pointsToMatric(x): longitude <-360

Does this mean that you have to use Geographic coordinate systems to use  distm?

Here is my code:

    proj4string(S1)

  [1] "+proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0
  +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

    head(coordinates(S1))

  coords.x1 coords.x2

  [1,]  330091.20  -3277120

  [2,] -672066.99  -3753558

  [3,] -671040.65  -3755774

  [4,] -676182.43  -3741369

  [5,]   61188.02  -3759233

  [6,]  156078.69  -3687183

    proj4string(S2)

  [1] "+proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0
  +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

    head(coordinates(S2))

  coords.x1  coords.x2

  [1,] -1246820.3  -586142.4

  [2,]  1831068.8  1024375.2

  [3,] -1093944.1  -617261.2

  [4,] -1075037.3  -614924.2

  [5,]  -652579.9 -3745182.5

  [6,]  -498183.2   261947.4

    test<-distm(coordinates(S1),coordinates(S2), fun=distHaversine)

  Error in .pointsToMatrix(x) : longitude < -360


texts:
What I want to achieve is to buffer a zip region geometry with a distance in meters and than get the coordinates of the new (expanded) polygon. The original geometry polygon is parsed from a kml file.

This is the code that I wrote so far:

Parse the kml file:

public Geometry getPolygon(String plz) {
  String poly = storage.get(plz);

  if (poly != null) {
    try {
      InputStream stream = new ByteArrayInputStream(poly.getBytes("UTF-8"));
      Parser parser = new Parser(new KMLConfiguration());
      SimpleFeature f = (SimpleFeature) parser.parse(stream);

      return (Geometry) f.getAttribute("GEOMETRY");
    } catch (UnsupportedEncodingException e) {
      e.printStackTrace();
    } catch (SAXException e) {
      e.printStackTrace();
    } catch (ParserConfigurationException e) {
      e.printStackTrace();
    } catch (IOException e) {
      e.printStackTrace();
    }
  }

  return null;
}

Transform the polygon:

public Geometry getBufferedZipPolygon(Geometry zipPolygon, double distanceMeters) {
  CoordinateReferenceSystem wgsCRS = null;
  CoordinateReferenceSystem distanceMetersCRS = null;
  try {
    wgsCRS = CRS.decode("EPSG:4326");
    distanceMetersCRS = CRS.decode("EPSG:32630");

    Geometry targetGeometry = zipPolygon.transform(distanceMetersCRS);

    Geometry buffer = targetGeometry.getBuffer(distanceMeters);
    buffer.transform(wgsCRS);

    return buffer;
  } catch (FactoryException e) {
    e.printStackTrace();
  } catch (TransformException e) {
    e.printStackTrace();
  }

  return null;
}

I have two Questions now. Is this transformation valid and how can I get the new Coordinates from my transformed Geometry Object?

Answer:

With the help from iant I finally figured out this solution.

Parser:

public Polygon getPolygon(String plz) {
  String polyString = storage.get(plz);

  if (polyString != null) {
    try {
      InputStream stream = new ByteArrayInputStream(polyString.getBytes("UTF-8"));
      Parser parser = new Parser(new KMLConfiguration());
      Polygon polygon = (Polygon) parser.parse(stream);
      return polygon;
    } catch (UnsupportedEncodingException e) {
      e.printStackTrace();
    } catch (SAXException e) {
      e.printStackTrace();
    } catch (ParserConfigurationException e) {
      e.printStackTrace();
    } catch (IOException e) {
      e.printStackTrace();
    }
  }

  return null;
}

Transformation:

public Geometry getBufferedZipPolygon(Polygon zipPolygon, double distanceMeters) {
  CoordinateReferenceSystem wgsCRS;
  CoordinateReferenceSystem distanceMetersCRS;
  try {
    wgsCRS = CRS.decode("EPSG:4326", true);
    distanceMetersCRS = CRS.decode("EPSG:32630", true);

    MathTransform transform = CRS.findMathTransform(wgsCRS, distanceMetersCRS);
    MathTransform reverseTransform = CRS.findMathTransform(distanceMetersCRS, wgsCRS);

    Geometry transformedZipPolygon = JTS.transform(zipPolygon, transform);
    Geometry bufferZipPolygon = transformedZipPolygon.buffer(distanceMeters);

    return JTS.transform(bufferZipPolygon, reverseTransform);
  } catch (FactoryException e) {
    e.printStackTrace();
  } catch (TransformException e) {
    e.printStackTrace();
  }

  return null;
}

Get the coordinates:

bufferedPlzPolygon = (Polygon) getBufferedZipPolygon(plzPolyStorage.getPolygon(zip), range);
Coordinate[] coords = bufferedPlzPolygon.getCoordinates();


texts: Reclassifying 16,000+ rasters and doing weighted sums on these rasters per year using ArcGIS ModelBuilder?

texts:
I've been trying to build a model in ModelBuilder (ArcGIS 10.1) which can do the following task:

I have a folder full of ca. 16,000 raster files (.asc) representing mean temperatures for each day of 40+ years with names like 19800101, 19800102 (first two days of January 1980) etc. Now I want:

import each .asc to raster
reclassify each raster (e.g. temperature x = 0, temperature y = 1)
do a weighted sum of all rasters of one year (to get the number of all days with temperature y)
write resulting yearly rasters to a different folder

I tried different things but nothing worked for the whole process. I can get 1) and 2) working in a model but fail to do 3) and 4) along with it automated and had to do it by hand which is very time consuming. I suppose it comes down to using different iterators, but I just can't find a way to get it working.

Can you point me in the right direction for doing this task automated and efficiently in ModelBuilder?


texts: Data defined label alignment option greyed out on new Quantum GIS Dufour

texts:
I am trying to customise the positioning of labels in relation to points in a vector file layer in the new QGIS 2.0.

In QGIS 1.8 I was able to enter two customised X/Y columns in the attribute table of a shapefile to precisely align the labels around the points.

However, from the new QGIS 2.0 menu 'Layer label options'>'Placement'>'Data defined', the X and Y 'Alignment' options are greyed out and I am unable to enter the columns.

Please see a screen dump of the issue below

Can anybody explain why this might be the case? Thanks for your consideration.


texts: How to dynamically rotate image in custom style in OpenLayers

texts:
I'm trying to make a custom Style in OpenLayers that replaces a path (Coordinate list) with an image at the end of the path pointing in the direction for the path (like the head of a snake). I am having trouble figuring out where/how to access the geometry to set the 'rotation' attribute in the Icon constructor. As geometry and renderer are the only options that allow for functions, I assume I will need to do my magic within a custom renderer. Are there default Icon render functions I could extend?

Here's my current code:

export class IconFrontPath extends Style {
  constructor(options: IFPOptions) {
    super(options);

    Style.call(this, {
      image: new Icon({
        src: options.url,
        rotation: ????
      }),
      geometry: this.geoFunction,
      zIndex: options.zIndex
    });
  }

  geoFunction: GeometryFunction = feature => {
    const coordinates = (feature.getGeometry() as LineString).getCoordinates();
    return new Point(coordinates[coordinates.length - 1]);
  };
}


texts: Rasterization of point vector in QGis: how to align the resulting raster to an existing one?

texts:
I have a point shapefile, with no associated .prj file, that I have to rasterize. Each point should be transformed into a pixel of 30x30 size aligned to the pixels of a Landsat image (the point file contains accuracy assessment points for a land cover semi-automated classification).

I don't know the original coordinate system but I'm reasonably sure of which one it is. The steps I took are:
1) tried saving the shp with different coordinate system, compared the results, and decided which might be the most likely one (because in the other cases the land cover classes assigned to the points would be mostly wrong).
2) saved the shp in that coordinate system, then converted into EPSG:32632 (same as the landsat images)
3) used gdal rasterize (from Qgis' Processing Toolbox) to rasterize the points. I set "output resolution in map units per pixel" as 30 x 30 m.

4) what I obtained are misaligned pixels like this:

The green pixel is the result of rasterization, the yellow point is the original point.

I tried redoing the whole procedure saving the original shp in different coordinate systems, but no luck.
I also tried the plugin rasmover to manually align the raster, but then it occurred to me that I don't know for sure how the raster should be aligned (i mean, I guess it should be aligned to the pixel it overlaps the most, but how can I be certain?). Moreover , if I use rasmover, many of the original points end up falling outside of the pixel that is supposed to contain them, because they are located near the edge.

So the questions are, what might be causing this problem? Is it related to the original georeferencing of the two files? And is there a way to fix this problem by using QGis, or by setting gdal_rasterize so that it will align the raster correctly? I've seen answers to similar problems related to ArcGis and GRASS but not specifically QGIS.


texts: Loop a Select by Attribute on all features (ArcPy)

texts:
I have a really basic script that selects an input Polygon by an Attribute and then does a series of Buffers on that, a select by location and then appends the result into a final output.  I need to run this process through each feature individually.  I created a tool parameter for the initial Select By Attribute, but how do I alter this script to loop through all features rather than manually altering the SQL statement?

Code below:

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# IndexMask_Standard_SQL.py
# Created on: 2015-09-02 09:01:40.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: IndexMask_Standard_SQL <Expression>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

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

# Local variables:
IndexPoly = "IndexPoly"
IndexPolys__3_ = IndexPoly
IndexPolys_Buffer5_Buffer_Mu1__3_ = IndexPolys__3_
IndexPolys_Buffer5 = "C:\\Users\\bwhiteley\\Documents\\ArcGIS\\Default.gdb\\IndexPolys_Buffer5"
IndexPolys_Buffer5_Buffer =  "C:\\Users\\bwhiteley\\Documents\\ArcGIS\\Default.gdb\\IndexPolys_Buffer5_Buffer"
IndexPolys_Buffer5_Buffer_Mu = "C:\\Users\\bwhiteley\\Documents\\ArcGIS\\Default.gdb\\IndexPolys_Buffer5_Buffer_Mu"
IndexPolys_Buffer5_Buffer_Mu1 = "IndexPolys_Buffer5_Buffer_Mu"
IndexMask_Export__3_ = IndexPolys_Buffer5_Buffer_Mu1__3_
IndexMask_Export = "C:\\Users\\bwhiteley\\Documents\\ArcGIS\\Default.gdb\\IndexMask_Export"

# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management(IndexPoly, "NEW_SELECTION", Expression)

# Process: Buffer
arcpy.Buffer_analysis(IndexPolys__3_, IndexPolys_Buffer5, "150 Feet", "OUTSIDE_ONLY", "ROUND", "NONE", "", "PLANAR")

# Process: Buffer (2)
arcpy.Buffer_analysis(IndexPolys_Buffer5, IndexPolys_Buffer5_Buffer, "5000 Feet", "OUTSIDE_ONLY", "ROUND", "NONE", "", "PLANAR")

# Process: Multipart To Singlepart
arcpy.MultipartToSinglepart_management(IndexPolys_Buffer5_Buffer,  IndexPolys_Buffer5_Buffer_Mu)

# Process: Make Feature Layer
arcpy.MakeFeatureLayer_management(IndexPolys_Buffer5_Buffer_Mu, IndexPolys_Buffer5_Buffer_Mu1, "", "", "OBJECTID OBJECTID VISIBLE NONE;SHAPE SHAPE VISIBLE NONE;PageNumber PageNumber VISIBLE NONE;ScaleRelative ScaleRelative VISIBLE NONE;ScaleDescription ScaleDescription VISIBLE NONE;Rotation Rotation VISIBLE NONE;PageName PageName VISIBLE NONE;SHAPE_Length SHAPE_Length VISIBLE NONE;SHAPE_Area SHAPE_Area VISIBLE NONE;BUFF_DIST BUFF_DIST VISIBLE NONE;ORIG_FID ORIG_FID VISIBLE NONE;SHAPE_length SHAPE_length VISIBLE NONE;SHAPE_area SHAPE_area VISIBLE NONE")

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(IndexPolys_Buffer5_Buffer_Mu1, "INTERSECT", IndexPolys__3_, "", "NEW_SELECTION", "INVERT")

# Process: Append
arcpy.Append_management("IndexPolys_Buffer5_Buffer_Mu", IndexMask_Export, "NO_TEST", "PageName \"PageName\" true false false 25 Text 0 0 ,First,#,IndexPolys_Buffer5_Buffer_Mu,PageName,-1,-1;SHAPE_Length \"SHAPE_Length\" false true true 8 Double 0 0 ,First,#,IndexPolys_Buffer5_Buffer_Mu,SHAPE_Length,-1,-1;SHAPE_Area \"SHAPE_Area\" false true true 8 Double 0 0 ,First,#,IndexPolys_Buffer5_Buffer_Mu,SHAPE_Area,-1,-1", "")


texts: Breaking single large raster into NumPyArrays using arcpy memory error

texts:
I'm trying to run raster calculator on a large raster file using the following approach : http://pro.arcgis.com/en/pro-app/arcpy/functions/rastertonumpyarray-function.htm, but I seem to run into memory errors.

I've used psutil and sys.getsizeof() to figure out that Python grows in memory usage each iteration at the myData line even though this object does not accumulate memory. It appears that Python is not releasing this memory and thus eventually maxes out its allocated memory.

I could research increasing memory allocation to Python, but think that there's got to be an easy way to get Python to release that precious memory. Here's the error code, but because I get it when Python reaches 1.4 GB RAM in task manager I think it's a memory issue.
 RuntimeError: ERROR 010240: Could not save raster dataset to blockraster\blockprocessingrdb22_191.tif with output format TIFF.

import arcpy
from arcpy.sa import *
#check out license
arcpy.CheckOutExtension("Spatial")
import os

# Loop over data blocks and reclass large raster
lake_elev = arcpy.sa.Raster(<LARGE RASTER>)
fileout = os.path.join(os.getcwd(),r"blockraster\blockprocessingrdb22.tif")

#settings
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = lake_elev
arcpy.env.cellSize = lake_elev

filelist = []
blockno = 0
blocksize = 1024

for x in range(0, lake_elev.width, blocksize):
    for y in range(0, lake_elev.height, blocksize):

        # Lower left coordinate of block (in map units)
        mx = lake_elev.extent.XMin + x * lake_elev.meanCellWidth
        my = lake_elev.extent.YMin + y * lake_elev.meanCellHeight
        # Upper right coordinate of block (in cells)
        lx = min([x + blocksize, lake_elev.width])
        ly = min([y + blocksize, lake_elev.height])
        #   noting that (x, y) is the lower left coordinate (in cells)

        #memory intensive operation
        myData = arcpy.RasterToNumPyArray(lake_elev, arcpy.Point(mx, my), lx-x, ly-y, nodata_to_value=-999)

        # Convert data block back to raster
        myRasterBlock = arcpy.NumPyArrayToRaster(myData, arcpy.Point(mx, my),
                                                 lake_elev.meanCellWidth,
                                                 lake_elev.meanCellHeight)

        # PROCESS DATA BLOCK -----------------------------
        # convert elevations to liklihoods
        myRasterBlock = Con(myRasterBlock<176.87, 0.9, Con((myRasterBlock >=176.87) & (myRasterBlock< 178.58), 0.8, Con((myRasterBlock >=178.58) & (myRasterBlock < 182.37), 0.7, Con((myRasterBlock >=182.37) & (myRasterBlock < 212.11), 0.6, 0.4))))
        # Save on disk temporarily as 'filename_#.ext'
        filetemp = ('_%i.' % blockno).join(fileout.rsplit('.',1))
        myRasterBlock.save(filetemp)

        # Maintain a list of saved temporary files
        filelist.append(filetemp)
        blockno += 1

# Mosaic temporary files
arcpy.Mosaic_management(';'.join(filelist[1:]), filelist[0])
if arcpy.Exists(fileout):
    arcpy.Delete_management(fileout)
arcpy.Rename_management(filelist[0], fileout)

# Remove temporary files
for fileitem in filelist:
    if arcpy.Exists(fileitem):
        arcpy.Delete_management(fileitem)

# Release raster objects from memory
del myRasterBlock
del lake_elev


texts: Draw line perpendicular to other line in ArcGIS Pro core host application

texts:
I am trying to draw a line perpendicular to a line segment in an ArcGIS Pro core host application.  There used to be a construct point in the ArcMap SDK but it is not there any more.  This is being done in .Net C#.
Is there  a tool for this that Esri has?


texts:
Currently I only have a GCS (WGS84), but it only allows for decimal degrees. Is there a global projection or CRS that would allow me to use Meters as a measurement?

I am using this for a tool that doesn't accept a linear unit (Value and Unit), but only a double and bases it off the CRS.


texts: Get rid of raster grid cells that are less than 1

texts: I'm using QGIS to work with raster data. I have a raster dataset of a country where each grid cell is assigned a number 1-100. However, the raster is in the format of a rectangle based on the extent of the file, so all of the cells outside of the borders of the country are assigned a 0. I cannot figure out how to create a GeoTIFF that excludes these cells that have a value less than 1. Any ideas? Thanks!


texts: Limiting decimal places in data imported as .CSV and joined to shapefile using QGIS?

texts:
I am importing a csv into QGIS 2.4 with a csvt file that specifies for percent columns Real(2.1) to limit it to 1 decimal place. It imports fine. When I join it to a shapefile, it looks fine. But when I save the shapefile with a new name and bring it back in (or use it in another program), all the numbers have about 10 decimal places. The extra 0s don't show, but some wind up as xx.x99999999999.

Is there a way of limiting the number of decimal places in those columns to just 1?


texts: Calculating spread of invasive species using ArcGIS for Desktop?

texts:
I have ArcGIS Desktop Advanced and so far have completed the following;

Used Supervised classification on a satellite image to identify the invasive species with field verification. I have converted the raster to a vector polygon to get the area in hectares.
Used Supervised classification in a second satellite image in the same area with a different time stamp and again converted to polygon to obtain area in hectares.

I would now like to;

Determine the rate of spread for this current period between the 2 satellite images/polygons and;
Use this to predict future spread of the species.

Can anybody offer some suggestions of the workflow to use?


texts:
I'm just starting to work with leaflet and have a basic question about showing a geoJSON file. First, I found this excellent article with example code:

http://bl.ocks.org/lxbarth/4019660

However I have a more basic/beginner question, in that example the set the JSON file through this url:

http://polymaps.appspot.com/county/{z}/{x}/{y}.json

However, the file I converted (from shapefile to geojson using quantum gis) is just a .geojson extension. I'm clearly missing something here, perhaps a way to convert the geojson file into a tile format? Can someone elaborate as to how I'm supposed to do this or what my missing step is? Do I need to be hosting the geoJSON file on a server in order to do this?

If my file is called polygons.geojson, how would I go about using it in the above reference code? Thanks!


texts:
I'm looking for an algorithm which, when given a particular route on a map with attributes such as incline/distance/shape/etc, can find a route which is similar (in terms of the attributes) but starts at a different point or in a different region on the globe.

Obviously it will be impossible in almost all cases to find a perfect fit but I'm looking for a "best match" sort of system with a method of measuring the similarity as well ideally.

I have tried looking, but most of my queries come up with either Map Matching problems or route similarity for GPS points along the same path. I may not know the correct terminology! Is there a name for this problem? What algorithm I can use to solve this?


texts: Excel Table Objects get invisible and lose data after re-opening map

texts:
I had a case when I had to Insert many Excel tables in the Map Layout and it worked perfectly, exported perfectly, etc. But when I saved it, closed the map and came back to it later ALL Excel tables disappeared. Moreover, if you click where the table should be, the table gets selected, but there is nothing besides the outline, no data, no background. Almost have a feeling that it got converted from Excel object to Draw object, since I can put background to it.
Any ideas, what could be wrong with it?
PS: I am using 10.2.2 ArcGIS and 2010 Excel


texts: Delete SDE feature class features with ArcPy

texts:
I'm working on a script that will empty my SDE database of all data. How would I go about deleting features from feature classes using ArcPy?

I have an SDE connection created with a user that should have delete privileges. So far I've tried using the arcpy.DeleteFeatures_management from code and from the Delete Features toolbox in data management.

My DB has SDE feature datasets with SDE feature classes inside.

The Delete Features tool box errors out with:

  ERROR 000732: Input Features: "Dataset Database Connections\Connection
  to sde.sde\User.MyFeatureDataset\User.MyFeatureClass" does not exist
  or is not supported.

The documentation is not really helping me. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000036000000


texts: Importing dgn geometry into empty file geodatabase feature class?

texts:
I have an empty geodatabase in arcgis for certain domain application. It has several feature datasets each containing related feature classes.
for example : feature dataset - NATURE,  it has a feature class  - TREEs

The geodatabase is empty.

Now i have features (corresponding to geodatabase) in a DGN drawing. My task is to import and load these features which are in the CAD drawing (.dgn file) to the corresponding features class in the geodatabase.

Can these kind of import and loading of data in database possible with FME. If yes how?

What information is required? I do not have any attribute information for features in .dgn file.

screenshot of dgn file. ( unique color and name)


texts: PyQGIS Standalone Crashes python.exe when zooming out too far / Process finished with exit code -1073741819 (0xC0000005)

texts:
I have the following PyQGIS Standalone App which is working fine so far.
One little problem is that if I zoom out to the MaxExtent(by scrolling the mousewheel) the app crashes python.exe so that the windows error report is startet.
Has anyone else encountered this problem and knows how to prevent it?

The error message is:

Process finished with exit code -1073741819 (0xC0000005)

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4 import *
from PyQt4 import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

# import Tkinter
# from Tkinter import *
import sys
#from qgis.core import QgsApplication
#from qgis.gui import QgsMapCanvas
import os

qgishome = r"C:/OSGeo4W64/apps/qgis"

app = QgsApplication([], True)
QgsApplication.setPrefixPath(qgishome, True)
QgsApplication.initQgis()
if len(QgsProviderRegistry.instance().providerList()) == 0:
    raise RuntimeError('No data providers available.')

canvas = QgsMapCanvas()
canvas.setCanvasColor(Qt.white)
canvas.enableAntiAliasing(True)
canvas.show()

layer = QgsVectorLayer("C:/daten/gwdat_31468.shp", "Grundwasserdaten", "ogr")
reg = QgsMapLayerRegistry.instance()
reg.addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayerSet([QgsMapCanvasLayer(layer)])
#QgsApplication.exitQgis()

app.exec_()

app.exitQgis()

Here's the version of QGIS i used:


texts: Contour polygon created at different intervals than the desired interval

texts:
I have a DEM in WGS84 coordinates with a minimum elevation of 2083 meters and a maximum elevation of 3000 meters. I then used the Contour tool in ArcGIS 10.7 to create contour polygons at an interval of 491.8 meters. Now the tool is creating contour polygons at an interval of 376 m instead of 491.8 meters. So, if you look at the first polygon, the contour max is 2459 meters, where as it should be 2,574.8 meters.
How can I solve this problem?
So that the contour polygons are created at the desired intervals of 491.8 meters.


texts:
I want to post-process GPS files via Differential processing (DGPS/DGNSS mode)

using rtklib (http://www.rtklib.com/).

I have rover and base observations in .rinex format.

I took FINAL satellite orbit solutions from NOAA-CDDIS archive (https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/orbit_products.html) in the following format:

WWWW/igsWWWWD.sp3.Z

Using the utility rtkpost from rtklib I get the error "no nav data".

Is that a problem of the sp3 format? (sp3 is listed within the acceptable formats)

Did I download completely wrong files for the satellite orbit solutions?


texts: Converting multiband-image to attribute table using ArcGIS Desktop?

texts:
I have a multiband-image with 12 Bands (monthly summed rain values, float)

How can I convert the multiband image into an attribute table in ArcGIS / ArcGIS Pro?


texts: Export Raster Information to Excel for statistical analysis

texts:
I want to make some statistical analyses on a raster. For this I do not need the spatial information. I just need to have a list with every cell value that I can open in Excel.

I tried exporting my raster as an ASCI xyz but it was for one too long for excel to open all values in one table and also the values were different from the ones that are displayed in my GIS.

Has anyone got a solution for this? Also I am open to any other solutions that lead to me being able to do complex statistical analyses (not only mean median etc.).


texts: OGR2OGR convert linestring to multilinestring

texts:
I have a shapefile with 3D linestrings. I want to convert this to 2D multilinestrings

I tried:-

ogr2ogr -t_srs EPSG:4326  Road_wgs84_2d_mlt.shp Road_wgs84.shp -nlt PROMOTE_TO_MULTI -dim 2 -progress

But this:-

ogrinfo -so -al Road_wgs84_2d_mlt.shp

Gives me:-

OGRFeature(Road_wgs84_2d_mlt1):0
  ID (String) = D291A5DD-66FA-4764-87F9-80D9FB8A83F5
  DISTNAME (String) = (null)
  ROADNUMBER (String) = B9097
  CLASSIFICA (String) = B Road
  DRAWLEVEL (String) = 0
  OVERRIDE (String) = F
  FEATCODE (Integer) = 15743
  LINESTRING (-3.444786910479451 56.164060208744381   52.480434944853187,-3.446986978754905 56.164483742630381 52.48251925688237)

As you can see the -nlt flag has not promoted the linestring to multilinestring. What am I doing wrong?

I have also tried -nlt MULTILINESTRING with the same result


texts: Building fishnet using polygon feature class as center of fishnet?

texts:
I'm interested in creating a fishnet where the origin of the fishnet builds from the center of a polygon feature class and builds out around it rather than building the fishnet starting from the lower left corner and only only the Y-axis and X-axis as depicted below. Using, ArcGIS, arcpy, python etc... is doable.

I've tried a variety of approaches, attempting to change the origin, at least think I was able to, or adding some sort of value to "buffer" or "pad" the input extent and or xy values, to no avail. I haven't come across much that can offer me in assistance thus far.

Creating fishnet grid with specified rotation angle in ArcGIS Desktop? includes a response:

  However, if you don't want to do the math or need to set your origin at a different point than the CRS origin, there's a step to take prior to starting the fishnet tool to get your y-axis coordinate.

However, I'm not sure how to implement.

I've read through:

Create Fishnet
How Create Fishnet works


texts: Finding node on polyline immediately after given point using Turf

texts:
I have 2 polylines, A and B. My objective is to

Find closest point on B from a given node on A
Find tangent on that point

So I use function nearestPointOnLine to find point on B closest to a node on A.
I am not sure how to achieve 2nd objective. I would probably need the node on B which comes immediately after the point I found for 1st objective, so that I have 2 points to determine the bearing.
Does Turf offer an API which can give me that node?


texts: Shapefiles and OpenStreetMap data do not overlap in QGIS?

texts:
So I'm pretty new to GIS, and I'm working on a mapping project using data for Vancouver, Canada. I want to combine OpenStreetMap data along with data provided by the City of Vancouver. But the City data don't line up with the OSM data when I map them in QGIS.

Following the advice of an answer to Projected maps should but do not overlap in ArcMap?, I checked the two data sources against a third (a map of the world). The OSM data is right -- British Columbia is in the right spot -- but the City data puts Vancouver in the south of France.

My data sources are:

Bikeways and libraries from the City of Vancouver. (Both are SHP files.) These are the data that don't align with the others.
A map of British Columbia from Open Street Map. I put the data into a database using osm2pgsql.
A world map that's part of the Open Street Style project.

Projections:

The City data are projected in UTM zone10, NAD83 (CSRS). I've selected this in QGIS's Layer Properties screen for those layers. Projection string: +proj=utm +zone=5 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
The OSM maps are in Mercator projection. QGIS has many Mercator choices, and I wasn't sure which one to use, so I guessed "Popular Visualization CRS / Mercator". +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
For the world map, Quantum GIS loaded the projection from the .prj file. That's +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs

Projected maps should but do not overlap in ArcMap? talks about zone problems. That sounds like what I might be experienceing. So I tried changing the libraries from UTM zone 10N to other zones (15N and 60N) but nothing changed.

Does anyone have suggestions for what to try next?


texts:
Edit: Yes, using a DEM for a bigger area worked. However, now I'm getting crazy slope values, everything is around 90 degrees. It's a hillside vineyard, I don't think it's vertical.

I am trying to create a slope from a clipped DEM. The DEM is for a small plot of land, as you can see below. Each square is 30x30m.

The slope map is not registering values, just going from NaN to NaN.


texts: Unique Values from Multiple Fields using Arcpy and Numpy

texts:
I am using Python 2.7 and am trying to generate a list of unique values from all values within two fields (cropsum and cropsum2). I am not trying to find unique row combinations. I've tried using both methods from http://geospatialtraining.com/get-a-list-of-unique-attribute-values-using-arcpy/ but neither method is working for me using two fields. My preference is to use numpy as and I am trying to write a script to manipulate data from a Frequency tool dbf output.

Method 1: Search Cursor

in_table = 'frequencySoutheastSubset.dbf'
field_names = ['cropsum', 'cropsum2']
with arcpy.da.SearchCursor(in_table, [field_names]) as cursor:
    uniqueCropCodes = sorted({row[0] for row in cursor})
    print uniqueCropCodes

returns the following error:
TypeError: 'field_names' must be string or non empty sequence of strings

Method 2: Numpy

arr = arcpy.da.TableToNumPyArray('frequencyNorthSubset.dbf',['cropsum', 'cropsum2'])
print(arr)

uniqueCropCodes = numpy.unique(arr)
print uniqueCropCodes

Returns unique row combinations, which is not what I want.
The numpy.unique documentation example

>>> a = np.array([[1, 1], [2, 3]])
>>> np.unique(a)
array([1, 2, 3])

suggests that the code above should work.

I tried suggestions from this np.unique(arr[['cropsum', 'cropsum2']].values)
did not work either: AttributeError: 'numpy.ndarray' object has no attribute 'values'

My dbf in case it's helpful:


texts: Figuring out if points are the same in FME

texts:
So I am working on a project using spatial point data in FME.  Both files contain 500 points and the hope is they are all the same.  To figure out if any are different I was thinking of using PointOnPointOverlayer tool instead of a spatial realtor.

What are your thoughts guys on this or is their a better way.  Just for heads up both files contain same information but from different sources so id and other useful information would not be the same to use a simple feature merger and then take unmerged ones as the different ones.

The goal is to figure out what records are not in the other data set.


texts: Using a dictionary to update one attribute based on another

texts:
I need to update a large dataset based on information in an Excel table. I figured the easiest way to do this would be using a dictionary I can fill out from the Excel table and run it on a feature class in a FGDB. Here is my code:

import os, arcpy

sitepoly = r'CodeTest'
AreaID = "907"##Change to current site

uCur = arcpy.UpdateCursor(sitepoly)
hs = uCur.next()
missingCnt = 0

while hs:
    OID = str(hs.getValue("OBJECTID"))
    gsID = str(hs.getValue("SITEMARKER_ID"))
    ##Define Datatypes from Excel using a dictionary here:
    ##Replace these values for different sections
    size= {"EE1":"3x8",
           "S":"3x8",
           "EE3":"5x10",
           "EE4":"3x8"}
    if gsID.startswith(AreaID):
        sectionID = str(hs.getValue("SECTION"))
        sizeattribute = str(hs.getValue("SITE_SIZE"))
        try:
            if len(sizeattribute) < 3 and sectionID == size:
                hs.setValue("SITE_SIZE", size)
                print str(gsID) + " has had the site size updated."
            uCur.updateRow(hs)
        except:
            missingCnt = missingCnt + 1
    hs = uCur.next()
del uCur

The code will finish, but "SITE_SIZE" is not updating with anything in the dictionary. For clarification: Each AREA has multiple SECTIONS and each SECTION has multiple features which I'm trying to update. Any ideas?


texts: Attempting to get the distance between two Point objects. Getting AttributeError: 'GeoSeries' object has no attribute '_geom'

texts:
I am a beginner in GIS. Using the code below to get the distance between two Point objects from two separate GeoDataframe. Using the shapely distance function.

Can someone please explain to me what I am doing wrong?

import geopandas as gpd

# Read file using gpd.read_file()
transit_data = gpd.read_file(transit)
stations_data = gpd.read_file(stations)
address_pts = gpd.read_file(city_pts)

# Point #1
address_pt = address_pts.loc[address_pts.full_address == '20 W 34th St, New York, NY 10001', 'geometry']

# Point #2
stop_pt = stations_data.loc[stations_data.name == '34th St - Hudson Yards', 'geometry']

# Calculate the distance between the 2 points
dist = stop_pt.distance(address_pt )

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-40-b60223ead16e> in <module>
----> 1 dist = stop_pt.distance(pt)

~\anaconda3\envs\gis\lib\site-packages\shapely\geometry\base.py in distance(self, other)
    440     def distance(self, other):
    441         """Unitless distance to other geometry (float)"""
--> 442         return self.impl['distance'](self, other)
    443
    444     def hausdorff_distance(self, other):

~\anaconda3\envs\gis\lib\site-packages\shapely\topology.py in __call__(self, this, other)
     44     def __call__(self, this, other):
     45         self._validate(this)
---> 46         self._validate(other, stop_prepared=True)
     47         d = c_double()
     48         retval = self.fn(this._geom, other._geom, byref(d))

~\anaconda3\envs\gis\lib\site-packages\shapely\topology.py in _validate(self, ob, stop_prepared)
     15
     16     def _validate(self, ob, stop_prepared=False):
---> 17         if ob is None or ob._geom is None:
     18             raise ValueError("Null geometry supports no operations")
     19         if stop_prepared and hasattr(ob, 'prepared'):

~\anaconda3\envs\gis\lib\site-packages\pandas\core\generic.py in __getattr__(self, name)
   5272             if self._info_axis._can_hold_identifiers_and_holds_name(name):
   5273                 return self[name]
-> 5274             return object.__getattribute__(self, name)
   5275
   5276     def __setattr__(self, name: str, value) -> None:

AttributeError: 'GeoSeries' object has no attribute '_geom'


texts: Add a tool to dropdown menu with the Python Add-In Wizard

texts:
Is it possible to add a tool created with the Python Add-In wizard to a menu in a toolbar that is also created with the add-in wizard?

I haven't had much luck in trying to figure this out.  I've searched around the internet/StackExchange and have not came up with anything.  Below is an example of what I'm trying to do.

I'm tyring to move the tool in the picture into the Tools dropdown menu

Here are the options I have when I right click on a menu

Here are the options I have when i right click on the toolbar.

Is it possible to edit the config.xml file and move the tool into a menu?


texts:
I have a 1 band 32 bit signed file system raster, with 4 attributes in addition to the value and count attributes with about 47000 rows, that is 55gb.

I wanted to reduce the size of this raster without losing spatial information and attribute information.


texts:
i have downloaded qgis-1.8.0.tar.bz2 and I want to install it on my linux computer(12.04 LTS). How can I perform it please?

I do it because I am desperate with Qgis 2.0.1 Dufour on ubuntu 12.04, when I wan to create a mapfile (I have many errors with RT-mapserver).

Thanks!


texts: How to create Race Route with arrow directions and distance markers using ArcGIS for Desktop?

texts: Anyone have any suggestions regarding how I would approach creating a race route with arrow directions and distance markers in ArcGIS? I could do it manually but I thought I recalled seeing something about a plug in that aided in this type of work.


texts: Postgis floating point error? in st_split


texts:
I am looking to find a way to simplify multi geometry KML files into single geometry files. However, rather than include the extra shape into my single geometry, I would like to remove it.

At the moment, I open the files in a text editor and manually remove the parts I need removed.

In the example below, I am looking to remove the two parts highlighted in yellow.

I create these files using ArcGis browser: create drive time area, then extract data based on city limits, and then manipulate them in batches in QGIS. The city limits I use are what create these multi-geometries.


texts: Show the coordinates at the centre of a map in OSGB format?

texts:
I am using QGIS v3.12 and projection EPSG:27700 - British National Grid.

I am trying to add a label to my map layout which will give the grid reference of the centre of the map as a UK 8 digit grid reference ( e.g. "AB123456" ).

So far I have managed to get the x and y coordinates using the expression:

Grid reference at map centre: [%round(x (map_get (item_variables('Map1'), 'map_extent_center')),0)%] , [%round(y (map_get (item_variables('Map1'), 'map_extent_center')),0)%]

Is there any function available in QGIS that can translate these coordinates into the required format?

Sorry my coding skills are not that great!


texts: Basemap layer disappear on print. Only happens with max zoom

texts:
I am fairly new to arcgis, please bear with me here.

ISSUE:

I uses standard esri.dijit.Print for map printing.
The base map displays normally in my map pane at all zoom levels.
The printing function works as expected for all except one zoom level.
When the map is zoom-in to the maximum level, the print result can show graphics from every layer except the base map. If only the base map and a single feature layer is selected, the result will show an mostly empty map with those feature graphics on it.
No javascript error and task response success.
tried comparing a successful print request of other zoom level with the failed request, no meaningful difference except the map extent.

lib version: 3.10.

connected via Lan

Request:
http://xx.xx.xx.xx/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute

Part of Web_Map_as_JSON: (note that the "land" layer is actually the base map I was referring to in the question)

"mapOptions":{"showAttribution":true,"extent":{"xmin":829211.7728597418,"ymin":840643.6650277233,"xmax":829322.236622336,"ymax":840675.9442589484,"spatialReference":{"wkid":102140,"latestWkid":2326}},"spatialReference":{"wkid":102140,"latestWkid":2326}},"operationalLayers":[{"id":"land","title":"land","opacity":1,"minScale":640000,"maxScale":1000,"url":"http://xxxx/MapServer"},{"id":"basemap","title":"basemap","opacity":1,"minScale":0,"maxScale":0,"url":"http://xxxx/MapServer"}

User-Agent:Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/47.0.2526.106 Safari/537.36

Does anyone have any idea where to even start with the debugging?


texts: GDAL C/C++ API and filtered pyramid tile generation

texts:
I have been tasked with adding functionality to an existing C++ developed toolset whereby a source dataset (typically a single GeoTIFF) needs to be split into tiles at various resolutions ('zoom' levels). Obviously this necessitates the generation of pyramid tile sets. Being completely new to the world of GIS I am finding there is quite a learning curve, hence would be grateful of any guidance with my particular problem.

So far I have developed a routine that generates the pyramid, and this is very much based on the gdal2tiles Python script. However, this made use of the RasterIO API calls and consequently the resulting tiles don't look so good when scaled down - i.e. very blocky and akin to a 'nearest-neighbour' sampling method. My question, is there a method (exposed via the  C/C++ API) that allows me to generate these tiles with bilinear or cubic filtering automatically? Or do I have to interrogate the band data myself and write the appropriate filtering kernels?

Due to the restraints imposed, I am unable to make use of the Python scripts and have to realise this functionality purely through the C/C++ API.


texts: Calculate meter offset around decimal degree position

texts:
I have a position given in decimal degrees (x.xxxxxxxx and y.yyyyyyyy). I need to draw a rectangle around it. The center of the rectangle matches the position. The dimensions of the rectangle is given in meters and it has a rotation ranging from 0-360 degrees.

Question
How can I calculate the four corners of the rectangle and return the result as four decimal degree values? Like arrayOf<LatLon> getRectangle(LatLon position, int rectWidthCm, int rectLengthCm, double rectRotation).

Example

I have a position given in LatLon format with two double values: latitude and longitude. We will assume this location is precice.
Then I get the task to draw a rectangle based on this position in a Google Maps chart. The rectangle can have any dimentions but let's use these in this example: Width = 0.9 meter and Length = 1.2 meters. Any heading may also be given so lets use this heading: 45. 0 Is north and going clockwise round (east = 90, south = 180 and west = 270). When the rectangle is pointing north it has the length in the north/south direction. Finally, the rectangle center should be equal to the given position.

Note: The project setup is an Android application with Kotlin support and a google maps chart. I am interested in a modern approach to this problem. Regarding precision loss it should at most be within centimeters.


texts: Add-In to Switch Selection [ArcMap 10.4]

texts:
I'm creating a suite of add-in tools on a tool bar, and using this tool as a test. How I want it to work is thus:

I have a buffer which has already been cut up by a line feature class
I select the portions of the cut-up buffer that I want to keep (will usually be
1-2)
I then press my add-in button, which switches my selection to all other parts of
the buffer (I will then add in 'delete all' functionality).

I have gone through the process and created the add-in project, toolbar, button etc., and loaded it into ArcMap through the customize menu. The button appears as intended, no 'Missing' issues. However, when I attempt to use it as described above, I get an egg-timer for a brief moment, and then nothing appears to happen. There is no change in the selection after I press the button compared to before.

My code:

import arcpy
import pythonaddins

class InverseSelection(object):
    """Implementation for RAFTv2_addin.invbtn (Button)"""

    def __init__(self):
        self.enabled = True
        self.checked = False

    def onClick(self):
        buffcut = "U:/ArcGIS/Default.gdb/buffcut"
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        arcpy.SelectLayerByAttribute_management(buffcut, "NEW_SELECTION", "INVERT")

    arcpy.RefreshActiveView()

I have also tried "SWITCH-SELECTION", but that didn't work either.


texts: Passing other arguments to event function in Leaflet

texts:
I am learning Leaflet. I want to create a function that will request WMS feature info when the map is clicked. And I want to do this for different maps, each with their WMS layers.

I imagine I will have to create a function that takes wms_url and map_object as parameters in the function in map.on('click', function). But I don't know how to pass parameters to the function in .on(). It seems .on automatically passes event object as a parameter to the function, and I don't see how to add other parameters.


texts: Translating GeoTIFF to 16-bit PNG with GDAL

texts:
I have a GeoTIFF file which displays the following information when running gdalinfo:

Band 1 Block=10012x1 Type=Float32, ColorInterp=Gray
  Description = Layer_1
  Min=1638.991 Max=2577.040
  Minimum=1638.991, Maximum=2577.040, Mean=2063.003, StdDev=185.546
  NoData Value=-3.4028234663852886e+038
  Metadata:
    LAYER_TYPE=athematic
    STATISTICS_COVARIANCES=34427.3350382379
    STATISTICS_MAXIMUM=2577.0400390625
    STATISTICS_MEAN=2063.0031636173
    STATISTICS_MEDIAN=0
    STATISTICS_MINIMUM=1638.9908447266
    STATISTICS_MODE=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=185.54604560119

Based on this information, I understand that the minimum and maximum grayscale values are 1638.991 and 2577.040 respectively.

What I want to do is convert this GeoTIFF to a 16-bit PNG image which shows ~938 (2577-1638) different shades of gray. I'm attempting to achieve this with the following gdal_translate command:

"c:\Program Files\GDAL\gdal_translate.exe" -outsize 4096 4096 -of PNG -ot UInt16 -scale 1638.991 2577.040 0 65535 "GDAL_IMG_2_GeoTIFF.tif" "GDAL_IMG_2_GeoTIFF.png"

My understanding here is that -ot UInt16 allows me to create a 16-bit image and -scale 1638.991 2577.040 0 65535 re-maps the geoTIFF grayscale range to the full 16-bit range (0 - 65535).

This command does indeed produce what appears to be a 16-bit image however when viewing the image properties in irfanView, the Number of unique colors is showing as 256.

Full bit-depth info:

Is there any step that I have missed here for creating 16-bit PNG's which make use of > 256 shades of gray or perhaps a better method?


texts: Raster to vector conversion in Google Earth Engine

texts:
I would like to convert a raster image into vector. My image is referenced to a specific zone (vect).
I'm working with Hansen dataset and in particularly the layer 'lossyear'. Raster image is from a different polygon in a specific year (for example 2007) and I obtain it by the function ReduceResolution. The value of different pixel that constitute the image are float number. The code is:

var loss_2007_vt = loss_2007_filtered_res.reduceToVectors({
  reducer: ee.Reducer.mean(),
  geometry: vect,
  crs: collection.projection(),
  scale:500,
  geometryType:'polygon',
});

I received this error message: "FeatureCollection (Error)
Image.reduceToVectors: Need 1+1 bands for , image has 1".


texts:
I would like to add a table with some statistics in the print composer.

I want to show the rate of feature for a specific value. Let say I have a "Existing" field with 'Yes' or 'No' properties. I want to add a table with Y% Yes / X% No

I tried to use the expression dialog but with no success.

For info: I am on mac OSx and can't use SQL requests as there is an important but since 2 years with this.


texts:
I have found some ArcObjects code to get unique values on one field for a IUniqueValueRenderer but I can't conceptually wrap my head around how to do that with 3 fields.  Would this be .Net code or is there a property/method in the IUniqueValueRenderer?

My goal is to get a list of layer files with python (these layer files have symbology from a master symbology layer which generally has more unique value combinations than what the target layer file has..meaning that some of the values in the IUniqueValueRenderer needs to be removed with the RemoveValue) then call an this arcobjects exe that will, once passed the path to a layer, finds it's unique values (from 3 fields) then check the renderer and start to remove the ones that don't apply.


texts: Adding WFS layer from GeoServer and OpenStreetMap using OpenLayers 2?

texts:
In a OpenLayers HTML file, I'd like to request the default tasmania data of my GeoServer (runs on localhost) and add them as a WFS Layer. The second layer should be an OSM layer.
This is my code up to now:

        function init() {

            var map;
            map = new OpenLayers.Map('map');
            map.addControl(new OpenLayers.Control.LayerSwitcher({'div':OpenLayers.Util.getElement('layerswitcher')}));
            map.addLayer(new OpenLayers.Layer.OSM());

            var roads = new OpenLayers.Layer.Vector("WFS", {
                strategies: [new OpenLayers.Strategy.BBOX()],
                protocol: new OpenLayers.Protocol.WFS({
                    url:  "http://localhost:8080/geoserver/ows?",
                    featurePrefix: "topp", //<-- GeoServer workspace name
                    featureType: "tasmania_roads",//<-- layer name
                    featureNS: "http://www.openplans.org/topp", //<-- Edit Workspace Namespace URI
                    srsName: "EPSG:4326",
                    geometryName: "the_geom" //<-- geometry field
                })
            });
            map.addLayer(roads);

            map.zoomToExtent(roads.getDataExtent());
        }

Concerning the OSM Layer, there are no problems in the result. But the data of the WFS layer is not being displayed. When I zoom to the WFS data extent (last code row) the map centers at (0,0) in the point of origin.

Could it have to do with projection problems?


texts: Merging and joining two line feature classes while maintaining attributes from both using ArcGIS for Desktop?

texts:
I created an event layer using the MAKE ROUTE EVENT LAYER tool based on anticipated scheduled events. The road layer which is the base layer is intact.

What I would like to achieve is combine layer 1 (base roads) with layer 2 (scheduled events) into one Result layer as shown in the figure below while maintaining attributes from both layers.

I have tried using Join, data join, spatial join, and other aspects and I still can't get it to do it right. If ArcGIS managed to merge the layers, it usually has missing data but for the most part my efforts have been unsuccessful.

what is the "best of practice way" to solve my problem?


texts:
I'm using Desktop 10.2 and I tried ,with this script, to add ecw file to mxd's but the result is a black & white raster(stretched symbology) instead of colored raster(RGB composite symbology) :

import arcpy
import arcpy.mapping
import os
import sys
from arcpy import env
import string

env.workspace = r"C:\project"
for mxd in arcpy.ListFiles("*.mxd"):
  print mxd
  mapdoc = arcpy.mapping.MapDocument(r"C:\project\\" + mxd)
  df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[0]
  Layer1 = arcpy.mapping.Layer(r"C:\project\layers\F1292_KVISH6_25CM_.ecw")
  arcpy.mapping.AddLayer(df ,Layer1 )
  mapdoc.save()
del mxd # i added

this is the image i get:

i don't know how to changed the python script so i will get RGB colored image ?


texts:
Goal: To export a labeled geo-referenced point layer created in ArcMap 10.2 to KML and have Google Earth only display labels of certain points when zoomed out past a certain scale (to avoid clutter).

I am tasked with creating a graduated circle map to be displayed in Google Earth. Here are the steps I took:

Exported table with values to ArcMap
Created graduated circles shapefile based on values in a column. Bigger values take on bigger circles
Labeled the circles based on another column of the attribute table. Labled with Mapplex label engine
Exported entire map to KML using map to KML conversion tool in ArcMap
Opened the KML file in Google Earth. Having all labels display at once clutters the map with labels

Questions:

How do I make it so that the labels for big circles appear when zoomed out and labels for smaller circles appear when zoomed in closer than a certain scale?

I tried doing this in ArcMap by setting the labels for smaller circles to not display when zoomed past a certain scale, but this doesn't affect the resulting KML file.

How do I make labels appear only when I hover my mouse over the point of interest in Google Earth? This is another option for me to reduce label clutter.

I am not that familiar with KML and don't know where to start.


texts: Global Human Settlement Layer-Rasterio and Fiona

texts:
I've downloaded the Global Human Settlement Layer data from here. They are huge files, too big for ArcMap or QGIS to handle comfortably so for that reason I decided to mask them. I was only interested in Ireland so I made sure my coordinate system was the same for both the tif and the polygon (Ireland-admin area) before I begun.

I tried to the following code on one of the tifs that I'm interested in.

import fiona
import rasterio
from rasterio import mask

with fiona.open(r"D:\GIS\Global Human Settlement Layer\Built_Eire\IRL_adm0_WEB.shp", "r") as shapefile:
    geoms = [feature["geometry"] for feature in shapefile]

with rasterio.open(r"D:\GIS\Global Human Settlement Layer\Built_Eire\GHS_BUILT_LDS2000_GLOBE_R2016A_3857_38_v1_0\GHS_BUILT_LDS2000_GLOBE_R2016A_3857_38_v1_0_p1.tif") as src:
    out_image, out_transform = mask(src, geoms, crop=True)
    out_meta = src.meta.copy()

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open(r"D:\GIS\Global Human Settlement Layer\Built_Eire\GHS_BUILT_LDS2000_GLOBE_R2016A_3857_38_v1_0\GHS_BUILT_LDS2000_GLOBE_R2016A_3857_38_v1_0_p1_Masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

The original .tif is 358.5MB, as soon as I try and run this code it started creating a huge file that even after only 9% completion it is 50GB.

Can anybody point me in the right direction as to what I may be doing wrong? I presume the file I downloaded was compressed?


texts: Disabling some topology rules temporarily to speed up checking other rules using ArcMap?

texts:
I am using ArcMap 10.6.

I have a topology with a large number of rules and many feature classes.  I need to check only certain rules and feature classes and wondering the best way to do this without having to validate the whole topology across all the rules?

I know I can search for only certain rules in the error inspector window, but I want to filter the process even before that, by telling the topology to ignore some rules until i turn them back on.

Is that possible?  Is there a better way to approach doing this when I have a large cadaster with many records and want to just check it against specific situations?


texts: Accessing SAGA algorithms (majority filter in particular) in QGIS

texts:
I'm using QGIS 3.18.1 and have downloaded SAGA 7.8.2 (following Integrating SAGA 7.8.1 into QGIS 3.16 (OS Windows 10)).
However, the majority filter option remains greyed out. I have tried installing the SAGA nextgen provider plugin, but that also appears to have no effect. What can I do to make the tool available?
Cheers
Ben


texts:
I am having difficulties completing a python tool.

I am trying to make an ID generator for unique values per grid. This will be a concatenated field of the layer (feature code, XX) + column of the fishnet (XX) + row of the fishnet (XX) + the Nth added feature (XXX)= the ID or 'FACILITYID'(9 digits, XXXXXXXXX)

Picture of the tool:

Picture of ArcMap

I.E., For Manholes I am trying to have the Feature Code be 60, and this will be the 106th manhole added into this cell. I would like the script to put the value of 601519106 into the selected NULL.

I keep getting this error.

line 44, in main
    if row.idField in nullvalues:
AttributeError: 'list' object has no attribute 'idField'

Here is my code so far:

arcpy.env.workspace = r"R:\Staging.gdb"

#Setting Parameters for tool
inputLayer = arcpy.GetParameter(0) #Input layer
i = arcpy.GetParameter(1) # Feature Class Code
idField = 'FACILITYID'
nullvalues = ["", " ", None, 0]
grid = r"R:\Toolbar\Shapefiles\Grids.shp"
arcpy.MakeFeatureLayer_management(grid, 'grid_lyr')
grid_lyr = 'grid_lyr'
selected_lyr = 'selected_lyr'
field_names = arcpy.UpdateCursor(grid_lyr)
cursor = arcpy.UpdateCursor(inputLayer)

#Select by NULL attribute
arcpy.SelectLayerByAttribute_management(inputLayer, 'NEW_SELECTION', '"FACILITYID" is NULL')
arcpy.SelectLayerByLocation_management(inputLayer, 'intersect', arcpy.SelectLayerByLocation_management (grid_lyr, 'intersect', inputLayer))
arcpy.SaveToLayerFile_management(inputLayer, selected_lyr, "relative")
ID = arcpy.GetCount_management('selected_lyr.lyr')
arcpy.SelectLayerByAttribute_management(inputLayer, "CLEAR_SELECTION")

for field in field_names:
    C = str(field.getValue('COL'))
    R = str(field.getValue('ROW'))

COL, useless1 = C.split('.')
ROW, useless2 = R.split('.')

#Starting the tool
arcpy.AddMessage("Generating ID's...\n")

FC = str(i).zfill(2)
Col = str(COL).zfill(2)
Row = str(ROW).zfill(2)
Fno = str(ID).zfill(3)
Code = str(FC + Col + Row + Fno)

with arcpy.da.UpdateCursor(inputLayer, idField) as cursor:
    for row in cursor:
        if row.idField in nullvalues:
            row.FACILITYID = Code
            cursor.updateRow(row)
            row = cursor.next()

del row
arcpy.Delete_management('grid_lyr')
arcpy.Delete_management(selected_lyr)
arcpy.AddMessage("\n\nFinished...")
arcpy.AddMessage("Last ID created: " + Code)


texts: Radius of curvature for curves in road segments for entire road network

texts:
I have a road network shapefile for a US state and I am trying to derive the radius of curvature for road segments in this shapefile. I need to be able to perform this task within just a few functions for thousands of road segments.

The road segments vary between 0.2 miles and maybe 5 or so miles. Each segment has several curves, but I am only interested in curves small enough to force a vehicle to make a reasonable turn. This should reduce the number of curves of interest per segment to a max of 1 or 2, with many segments not having such curves at all. I have not determined the maximum radius of the curve (the threshold) I would like to use.

For segments with only 1 such curve, I would like to obtain the values of those curves for those segments, and for those with more than 1, an average will do.


texts:
I'm new to using GRASS GIS and GIS in general. I have a point layer and a line layer which I added on to the layer tree using two .shp files. I need to use v.net tool to create a network by connecting these two layers as edges and nodes respectively so that I can later extract the Steiner tree using v.net.steiner. So I typed in the following code into the console in order to do this:

v.net input=cmbRoads@steiner points=bsSnapped@steiner
out=network_net \ operation=connect threshold=500

All I get is the following error message:

I am using GRASS GIS 7.6.0 on Windows 10.


texts: Multiple Geoserver GetFeatureInfo templates for same layer?

texts: I have layer 'states' on Geoserver for which I need multiple GetFeatureInfo templates. I would like to use that layer a couple of times on different maps, but I want different data in Leaflet popup. Is it possible to define multiple templates for same layer? Is it possible to select template name in GetFeatureInfo request? Is there any way to select template?


texts: What is the correct way to convert a USGS Orthoimage from WorldImage format to tiled GeoTiff for GeoServer?

texts:
USGS provides different resolution orthoimage quadrangles in the WorldImage format (*.tif with an accompanying *.tfw world file). Some of these quadrangles can be used by GeoServer (2.1.0 - 2.1.2) by creating a WorldImage store and publishing the layer (some of them throw an exception when trying to publish the layer). When publishing a GeoTiff layer that has been converted from a currently working WorldImage layer, I am getting the following exception:

2012-01-10 08:38:33,228 WARN [referencing.factory] - Axis elements found in a wkt definition, the force longitude first axis order hint might not be respected:
PROJCS["WGS84 / Simple Mercator", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295]], PROJECTION["Mercator_1SP_Google"], PARAMETER["latitude_of_origin", 0.0], PARAMETER["central_meridian", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH], AUTHORITY["EPSG","54004"]]
2012-01-10 08:38:39,674 INFO [geoserver.web] - Error saving layer
java.lang.NullPointerException
    at org.geoserver.catalog.impl.CoverageStoreInfoImpl.getFormat(CoverageStoreInfoImpl.java:42)
...

I tried using gdal_translate in the following way (substitute 1234 with real values): gdal_translate -of GTiff -a_srs "EPSG:2239" -co "TILED=YES" 1234.tif ../geotiff/1234.tif The conversion succeeds, and I can successfully create a GeoTiff store in GeoServer, but cannot publish the layer because of the error described above. This leads me to believe that I am trying to publish layers (in both cases) that are not properly formatted. Thus, what is the correct way to convert a [working] USGS orthoimage in WorldImage format to a tiled GeoTiff for GeoServer 2.1.2?


texts:
I'm using QField and was able to setup almost everything.
I managed to enable photos creating a table (with the "attachment" type, relative path and so forth) and building relations to my Shapefile. This all works great.

My problem is that when I take a photo from the shape's menu, I get weird behavior:

Sometimes, the photo is stored in DCIM at root of the project in the
phone
Sometimes, the photo is stored in the DCIM of the file, outside the project folder and mixed with my regular phone photos
In that second situation, the photo's name is different from the name stored in my table

Is there a way to force where the pictures taken from QField are stored?
I'm not using a external memory card, everything is stored locally on the phone.

Here are some pictures of my project:


texts: QGIS 3 modeler: how to use the fields of an intermediate output as input in another algorithm?

texts:
I'm working on a model to calculate land uses (usos) in different areas (zonas):

I have corrected geometries, cutted the uses layer, corrected geometries again and used the fields calculator to assing the area ID to every uses polygon contained in it [geomintersects ('zonas', 'ID')]. After that, I want to use the "collect geometries" algorithm to summarize uses by code and area. However, none of the fields of this last output appears in the dialog box.

I have searched through the "input" window of the modeler without any success. How can it be done?


texts: ArcGIS for WPF: LayerCanvas.ArrangeOverride Setting Invalid FrameworkElement Size

texts:
I'm reposting this here from ESRI's ArcGIS for WPF Forum in hopes of getting some more exposure. Original post is here.

I'm intermittently experiencing an error when rapidly changing map scales on an ArcGIS WPF Map Control (v. 10.1.1) (i.e. from very large to very small, or vice-versa) wherein my application will crash. It appears that the LayerCanvas.ArrangeOverride is, for some reason, attempting to set an invalid Width property on the FrameworkElement's dependency property. Based on this link, it appears that the issue was carried over from the Silverlight API.

Has anyone else encountered this problem and been able to work around it?

The following is a snippet of the stack trace with the offending method:

System.ArgumentException: '-907.750834278924' is not a valid value for property 'Width'.
   at System.Windows.DependencyObject.SetValueCommon(DependencyProperty dp, Object value, PropertyMetadata metadata, Boolean coerceWithDeferredReference, Boolean coerceWithCurrentValue, OperationType operationType, Boolean isInternal)
   at System.Windows.DependencyObject.SetValue(DependencyProperty dp, Object value)
   at System.Windows.FrameworkElement.set_Width(Double value)
   at ESRI.ArcGIS.Client.LayerCanvas.ArrangeOverride(Size finalSize)
   at System.Windows.FrameworkElement.ArrangeCore(Rect finalRect)
   at System.Windows.UIElement.Arrange(Rect finalRect)
   at System.Windows.ContextLayoutManager.UpdateLayout()
   at System.Windows.ContextLayoutManager.UpdateLayoutCallback(Object arg)
   at System.Windows.Media.MediaContext.InvokeOnRenderCallback.DoWork()
   at System.Windows.Media.MediaContext.FireInvokeOnRenderCallbacks()
   at System.Windows.Media.MediaContext.RenderMessageHandlerCore(Object resizedCompositionTarget)
   at System.Windows.Media.MediaContext.RenderMessageHandler(Object resizedCompositionTarget)
   at System.Windows.Threading.ExceptionWrapper.InternalRealCall(Delegate callback, Object args, Int32 numArgs)


texts:
I'm rendering maps into SVG files, and I'd like to put geographic metadata into those files, as described in section 7.13 and 7.14 of the SVG spec, but I'm uncertain about how to interpret the <svg:transform> element.

My primary use case for this is to eventually get the information into a GeoTIFF file. Since I have full control over the rendering chain, I could convey this information by any sneaky means that I care to invent, but since there is a documented way of doing it...

Let's take a concrete example:

I'm drawing my map on A4-sized paper, and I have set up the viewBox so that the viewBox coordinate system gives me 90 units per inch. In other words, the top left corner of the paper is at (0, 0) and the lower right corner is at (744, 1052).

My geographic data is in EPSG:3006 (SWEREF99TM, more or less equivalent to EPSG:32633), so a typical point P could be (660000, 6612000). Possibly, those coordinates should be flipped, but that's a separate can of worms -- this is what works with PostGIS, so it is what I have.

I want a scale of 1:50000, and the lower left corner of the map at (655000, 6606000). This puts my point P roughly at the centre of the page.

I want a sane coordinate system on my page, so I wrap the entire map in

<g transform="translate(0 1052) scale(1 -1)">

thus putting the origin at lower left, and Y axis upwards.

I then compute coordinates for my point in this coordinate system:

x = ( 660000 -  655000) / 50000 * 90 / 0.0254 = 354
y = (6612000 - 6606000) / 50000 * 90 / 0.0254 = 425

And finally I draw my point as a small circle

<circle cx="354" cy="425" radius="3" />

(In reality, I keep a few more decimals.)

Now, what should I put in <svg:transform> to describe this? Should I describe the transform from geographic coordinates to the actual numbers written into the file, i.e.

scale(0.0709) translate(-655000 -6606000)

or should I describe the transform to viewBox space

translate(0 1052) scale(1 -1) scale(0.0709) translate(-655000 -6606000)

or something completely different?


texts: Drawing ellipses over a tile map server in PyQGIS?

texts:
I have a few ellipses, defined by a latitude, longitude, semi-major axis in nautical miles, semi-minor axis in nautical miles, and an orientation. I'm trying to render them over a Google Maps layer.

I create the QgsEllipse:

def qgs_ellipse(self):
    return QgsEllipse(QgsPoint(self.latitude, self.longitude), self.semi_major_axis, self.semi_minor_axis, self.orientation)

Create the QgsFeature:

def qgs_feature(self):
    feature = QgsFeature()
    feature.setGeometry(self.qgs_ellipse().toPolygon())
    return feature

Create a QgsRasterLayer from a Google Maps TMS:

self.tms_layer = QgsRasterLayer("type=xyz&url=" + tms[1], tms[0], "wms")
self.project.addMapLayer(self.tms_layer)

Create the QgsVectorLayer:

self.ellipse_layer = QgsVectorLayer("polygon", "Ellipses", "memory")
self.ellipse_provider = self.ellipse_layer.dataProvider()
self.project.addMapLayer(self.ellipse_layer)

Add the ellipses to the vector layer:

def add_ellipses(self, ellipses):
    features = []
    for ellipse in ellipses:
        features.append(ellipse.qgs_feature())
    self.ellipse_provider.addFeatures(features)

And finally add the layers to the QgsMapCanvas:

self.canvas.setLayers([self.ellipse_layer, self.tms_layer])
self.canvas.setExtent(self.ellipse_layer.extent())
self.canvas.show()

The issue is that all the ellipses show up in the center of the map, all centered on the same point. I also don't know if I need to convert the ellipse axis' from nautical miles to something else.

Is there a coordinate system or something I need to set up?

I am pretty new to GIS.


texts: Leaflet TypeError: this._northEast is undefined

texts:
I have some codes used Leaflet and I want to refactor them below one is working:

jQuery(document).ready(function(){
                    $.post("StatisticsResultServlet",
                    {
                      queryId:queryIdStr,
                      field:fieldName,
                      groupingLevelValue:groupingLevel,
                      selectedCityName:selectedCity,
                      divideByPopulation:divideByPop,
                      nofLegendGap:numberOfLegendGap,
                      multipliedNormalizationConstant:multipliedConstant
                    },
                    function(data,status){

                        statisticsResultArray = data.resultArray;
                        legendRangeArray = data.legendRanges;

                        var geojsonFeature;
                        if(groupingLevel == 1){
                            geojsonFeature = BOLGE;
                        }else if(groupingLevel ==2){
                            geojsonFeature = IL;
                        }else if(groupingLevel ==3){
                            geojsonFeature = eval(selectedCity);
                        }
                        geojson = L.geoJson(geojsonFeature, {
                        style: style,
                        onEachFeature: onEachFeature
                        });

                        geojson.addTo(map);

                        var southWest  = geojson.getBounds().getSouthWest();
                        var northEast  = geojson.getBounds().getNorthEast();
                        var bounds = new L.LatLngBounds(southWest, northEast);
                        map.fitBounds(bounds);

                        prepareLegend();

                    });

                });

This method is too long and I separated to different classes the business logic. Related my class which I got the error is below:

function PolygonService(){
     this.createPolygon = function(geoJsonFeature){....
     this.geoJsonBounds = function(map){
        var boundsOfGeoJson;

        var northEast  = geoJson.getBounds().getNorthEast();
        var southWest  = geoJson.getBounds().getSouthWest();
        boundsOfGeoJson = new L.latLngBounds(southWest, northEast);

        map.fitBounds(boundsOfGeoJson);

    };

The error I got

  TypeError: this._northEast is undefined

How to solve the problem?


texts:
I want to extend metadata style. I want to add three new fields in metadata of feature class.

I know there are different standards of metadata which ArcGIS Supports.

Is there any way or possibility to have my own standard or customize some existing and made it available in ArcGIS in my organization?

I know that there is a way to create an xml file in catalog and view it as meta data and do the changes and use it as template, but can i add new fields even to this template?

If yes, Can i use that in ArcGIS Portal too ?


texts: Why doesn't my Mapbox + Leaflet viz look retina-friendly?

texts:
This is probably a very simple question but it's not clear to me.

The Mapbox docs say:

  When you create a project from the editor, Retina (high resolution)
  tiles are automatically used on devices that supports them. If you are
  using Mapbox.js, retina tiles are always automatically used when
  available. No extra code is necessary.

Yet my projection looks fuzzy compared with the map in my Mapbox Projects section.

I use mapbox.js and .css:

<link href='https://api.tiles.mapbox.com/mapbox.js/v1.0.0/mapbox.css' rel='stylesheet' />
<script src='https://api.tiles.mapbox.com/mapbox.js/v1.0.0/mapbox.js'></script>

Here are pictures of the respective maps:

The projection (made with Leaflet):

The project as it appears on Mapbox.com's "Project" page:

What mistake am I making?


texts: import SpatialPolygonsDataFrame (SPDF) into PostgreSQL DB schema from R

texts:
This is a followup question to
Writeogr to PostgreSQL/PostGIS database with R

How can I choose the schema where to put the SPDF (or any other spatial object)?

postgis_insert(con, df=catchments, tbl="land_use", geom_name="st_astext")

I can't just type schema1.land_use in the tbl argument and there is no further argument to assign the schema. It seems that I can only put it into the public schema of the DB (which is kind of the default schema in PostgreSQL)


texts: Exporting Transparent Background from QGIS?

texts:
I have QGIS version 3.4.3 which does not have print composer.  I have a series of lines/tracks/boundaries, etc that are in a layer.  I am trying to export that layer with a transparent background for installation in a Garmin GPS unit.  I have successfully completed all steps with the exception of creating the transparent background.

I have scoured the web, looked through the technique in Print Composer which is not available in this version and tried to no avail at GPS File Depot.

What am I overlooking?